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Abstract 


A  system  of  equations  governing  the  incompressible  flow  through  a  multiple-blade-row 
turbomachine  is  presented.  These  equations  basically  form  the  pseudo-compressible  analog  to  the 
equations  derived  by  Adamczyk  [1984]  for  compressible,  average-passage  flow.  The  methodology 
behind  the  derivation  is  outlined,  including  a  closure  model  for  the  time-averaged  form  of  the 
equations.  The  equations  are  then  preconditioned  to  facilitate  numerical  treatment.  An  explicit 
numerical  procedure  based  on  Runge-Kutta  time  stepping  for  cell-centered,  hexahedral  finite 
volumes  is  outlined  for  the  approximate  solution  of  the  governing  equations.  Convergence 
acceleration  techniques,  boundary  conditions,  and  closure  issues  are  also  addressed  for  the 
numerical  scheme.  Finally,  results  are  presented  for  a  simulation  of  the  high  Reynolds  number 
flow  through  a  two-blade-row,  axial-flow  pump.  These  comparisons  suggest  that  the  pseudo- 
compressible  average-passage  equations  can  make  reasonable  predictions  of  the  highly  three- 
dimensional  flow  within  a  multiple-blade-row  turbomachine  operating  in  an  incompressible  flow 
regime.  However,  especially  in  wake  regions,  it  is  clear  that  the  behavior  of  the  algebraic  eddy 
viscosity  model--at  least  with  the  present  grid-requires  improvement  for  the  accurate  prediction  of 
the  evolution  of  the  downstream  velocity  field. 
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1.  Introduction 


In  recent  years,  the  design  and  analysis  of  turbomachinery  has  increasingly  turned  to 
numerical  computations  in  order  to  improve  estimates  of  this  complex  flow  field.  The  improved 
estimates  of  this  three-dimensional,  viscous  (turbulent)  flow  becomes  even  more  important  when 
dealing  with  turbomachines  with  a  multiple  number  of  blade  rows.  Most  of  the  recent  advances  in 
this  area  have  dealt  with  compressible  flow.  This  report  summarizes  work  to  extend  numerical 
computations  of  multiple-blade-row  turbomachines  to  incompressible  flow. 

First,  this  report  outlines  the  mathematical  basis  of  the  average-passage  representation  of  a 
multiple-blade-row  turbomachine.  Specifically,  starting  with  the  Navier-Stokes  equations, 

Chapter  2  describes  the  basic  methodology  that  leads  to  the  incompressible  average-passage 
equations.  Then,  the  closure  problem  associated  with  this  representation  is  discussed,  followed  by 
a  description  of  a  closure  model  appropriate  for  a  single-stage  (two-blade-row)  turbomachine 
operating  incompressibly.  Finally,  the  governing  equations  are  cast  in  a  preconditioned,  pseudo- 
compressible  form  that  is  more  appropriate  for  numerical  approximation. 

After  establishing  the  mathematical  model,  we  describe  the  numerical  technique  used  to 
approximate  a  steady-state  solution  to  the  system  of  equations  presented  in  Chapter  2.  Specifically, 
Chapter  3  describes  the  spatial  and  temporal  discretizations  of  the  governing  equations,  along  with 
the  acceleration  techniques  used  to  expedite  the  achievement  of  a  steady-state  solution  to  the 
discrete  equations.  Then,  we  give  a  brief  description  of  the  boundary  conditions  appropriate  for  a 
pseudo-unsteady,  average-passage  representation  of  a  blade  row.  Finally,  Chapter  3  addresses  the 
numerics  of  the  closure  model  proposed  in  Chapter  2  for  a  single-stage  machine. 

With  the  development  of  a  computer  code  called  ISTAGE  that  numerically  solves  the 
average-passage  equations  for  an  incompressible  flow,  we  computed  an  experimentally-measured 
flow  field  as  a  first  step  towards  code  validation.  Zierke,  Straka,  and  Taylor  [1993]  acquired 
measurements  in  the  high  Reynolds  number  pump  (HIREP)  facility,  during  an  experimental 
program  that  was  performed  in  parallel  with  the  numerical  program  described  herein.  The 
numerical  computations  of  this  incompressible  flow  of  water  through  a  two-blade-row  pump  were 
performed  as  a  prediction,  using  a  uniform  inflow  rather  than  the  experimentally  measured  inflow. 
Chapter  4  shows  a  large  number  of  comparisons  between  the  numerical  predictions  and  the 
experimental  data,  resulting  in  a  good  measure  of  the  current  status  of  our  numerical  model. 
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2.  Governing  Equations 


In  this  chapter,  the  derivation  of  the  average-passage  equation  system  is  briefly  described. 
In  his  original  derivation,  Adamczyk  [1984]  outlined  a  more  rigorous  treatment  for  compressible 
flow.  Beginning  with  the  Navier-Stokes  equations,  we  describe  the  general  form  of  the  average- 
passage  equations  for  incompressible  flow,  including  the  averaging  operators  that  are  used  in  their 
derivation.  This  is  followed  by  a  description  of  a  simple  closure  model  to  estimate  the  resulting 
body  forces  and  correlation  terms.  Finally,  the  equations  are  cast  in  a  preconditioned  pseudo- 
compressible  form  better  suited  for  numerical  treatment. 


2.1  Navier-Stokes  Equations 


The  starting  point  for  the  derivation  of  the  average-passage  equations  is  the  incompressible 
Navier-Stokes  equations  in  cylindrical  coordinates.  They  are  given  as  follows: 
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Respectively,  these  equations  represent  the  conservation  of  mass,  axial  momentum,  radial 
momentum,  and  angular  momentum.  The  thermal  energy  equation  is  not  included  because  it  is 
decoupled  from  the  above  equations  if  the  kinematic  viscosity  is  assumed  to  be  constant. 
Equations  2.1  through  2.4  are  written  in  terms  of  nondimensionalized  dependent  and  independent 
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variables  as  discussed  later.  The  shear  stress  terms  in  the  Navier-Stokes  equations  are  given  by 
the  following: 
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where  the  reference  Reynolds  number  is  defined  as 
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For  the  above  equations,  the  following  nondimensionalizations  apply: 
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with  the  assumption  of  incompressibility  giving 


P 


P 

P  ref 


1.0 


and 


1.0  . 


For  completeness,  the  shaft  rotation  rate  is  nondimensionalized  as 


fi  =  reS 


U. 


rtf 


In  the  above,  a  tilde  (~)  denotes  a  dimensional  quantity  and  the  subscript  "ref'  denotes  a  reference 
quantity.  In  this  work,  the  equations  are  applied  to  the  flow  through  a  turbomachine.  The 
reference  length  is  taken  to  be  the  machine  diameter  and  the  reference  velocity  is  taken  to  be  the 
average  velocity  through  the  inlet  of  the  machine  (that  is,  the  mass  flow  divided  by  the  inlet  area). 
The  reference  density  and  viscosity  are  taken  from  the  properties  of  the  working  fluid. 

Equations  2.1  through  2.4  represent  the  starting  point  for  all  of  the  analysis  that  follows  in 
the  remainder  of  this  report.  Note  that  in  their  present  form,  these  equations  are  exact— they  are 
the  equations  governing  the  unsteady  motion  of  an  incompressible  fluid.  However,  because  of  the 
lack  of  a  time  derivative  in  Equation  2.1,  the  equations  cannot  be  written  as  a  complete  system  of 
time  evolution  equations  (unlike  the  equations  governing  compressible  flow).  Ultimately,  this 
difficulty  will  be  rectified  by  restricting  the  interest  to  steady-state  solutions,  permitting  the  familiar 
pseudo-compressibility  assumption.  However,  this  is  done  only  after  the  derivation  of  the  average- 
passage  equations  of  motion. 


2.2  Average-Passage  Equations 

The  average-passage  equations  govern  the  "average"  flow  within  a  blade  passage  embedded 
in  a  multiple-blade-row  turbomachine.  Deriving  equations  that  govern  an  averaged  flow  field  is  a 
familiar  exercise.  Ensemble  (or  Reynolds)  averaging  the  Navier-Stokes  equations  to  yield  the  so- 
called  Reynolds-averaged  Navier-Stokes  (RANS)  equations  is  a  well  known  procedure  used  to 
formulate  equations  describing  turbulent  (nondeterministically  unsteady)  flows  in  some  "average" 
sense.  In  that  case,  the  term  "average"  means  steady  or,  at  most,  deterministically  unsteady.  For 
the  average-passage  equations,  the  term  "average"  is  understood  to  be  flow  that  is  steady  and 
spatially  periodic  over  the  pitch  of  the  blade  row  of  interest.  Clearly,  if  a  turbomachine  has  M 
blade  rows,  there  are,  in  general,  M  associated  average-passage  flow  fields;  that  is,  each  blade  row 
has  a  representative  or  average  passage. 

Adamczyk  [1984]  developed  the  average-passage  equations  by  the  sequential  application  of 
three  distinct  averaging  operators  to  the  Navier-Stokes  equations  within  a  multiple-blade-row 
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environment.  The  resulting  equations  thus  describe  the  flow  through  a  referenced  blade  row 
embedded  in  a  multiple-blade-row  turbomachine  in  a  steady,  spatially  periodic  sense  that  is  unique 
to  the  referenced  blade  row.  A  difficulty  associated  with  this  averaged  system  of  equations, 
however,  is  the  necessity  to  provide  body  forces  and  temporal  and  spatial  correlations  to  properly 
account  for  the  effects  of  neighboring  rotating  and  stationary  blade  rows  on  the  referenced  blade 
row.  Analogous  to  Reynolds  averaging,  the  correlations  are  a  direct  result  of  the  nonlinear 
convection  terms  in  the  governing  equations  and  must  be  modelled  by  additional  equations  or 
empirically-based  approximations.  These  operators  are  briefly  described  in  this  section. 


2.2.1  Ensemble  Averaging 

The  first  operator  applied  to  the  governing  equations  (2.1  through  2.4)  is  the  familiar 
ensemble  (or  Reynolds)  averaging  operator.  For  an  arbitrary  function,  / ,  this  averaging  operation 
is  defined  as 


f  =  lim  Tj  E  fi 

N-*oo  N  /  =  1 


(2.5) 


where  fi  is  the  ith  realization  of  / .  The  function  can  then  be  decomposed  into  an  ensemble- 
averaged  component  plus  a  random  component,  that  is 

/-/♦/  (2'6) 


where,  by  construction, 


It  is  important  to  note  that  the  ensemble  averaging  operator  commutes  with  both  temporal  and 
spatial  differentiation. 

The  procedure  for  ensemble  averaging  the  governing  equations  begins  with  decomposing 
the  velocities  and  pressure  in  Equations  2.1  through  2.4  according  to  Equation  2.6.  The  equation 
is  then  operated  on  by  Equation  2.5,  taking  advantage  of  its  commutative  properties  with 
differentiation.  The  result  is  the  familiar  Reynolds-averaged  Navier-Stokes  equations.  The  RANS 
equations  are  merely  the  Navier-Stokes  equations  with  the  dependent  variables  replaced  by  their 
ensemble-averaged  counterparts  and  with  the  inclusion  of  some  additional  stress-like  terms.  These 
apparent  or  "Reynolds"  stresses  are  a  direct  result  of  the  nonlinear  convection  terms  and  involve 
correlations  of  the  randomly  fluctuating  components  of  velocity.  To  solve  the  RANS  equations 
requires  some  approximation  or  modelling  of  these  terms.  This  is  the  familiar  closure  problem  of 
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turbulent  flow  prediction.  There  are  in  existence  a  myriad  of  models  available  for  Reynolds  stress 
closure.  A  fairly  standard  model  is  utilized  in  this  work  and  is  described  in  the  next  chapter. 

In  essence,  by  solving  the  RANS  equations  (after  providing  some  appropriate  model  for  the 
Reynolds  stresses),  a  deterministic  flow  field  is  sought  with  the  effects  of  the  random  fluctuations 
in  the  variables  accounted  for  in  some  average  sense  Conceptually,  this  is  what  is  intended  for 
two  additional  averaging  operations:  time  averaging  and  passage-to-passage  averaging. 

2.2.2  Time  Averaging 

The  second  step  in  the  derivation  of  the  average-passage  equations  is  the  application  of  a 
time-averaging  operator  to  the  ensemble-averaged  governing  equations.  The  operator  is  designed 
so  that  it  may  be  applied  everywhere  in  the  flow  domain  including  the  regions  within  rotating  blade 
rows.  Basically,  this  operator  averages  its  operand  over  a  period  equal  to  one  shaft  rotation.  The 

form  of  this  operator  acting  on  an  ensemble-averaged  function,  /  ,  is  given  by 


(.  *2t/Q 


(2.7) 


where 

Q  =  shaft  speed 

XR  =  geometric  blockage  of  the  neighboring  rotating  blade  row(s) 

H(t)  =  "gate"  function. 

The  purpose  of  the  gate  function  is  to  permit  the  application  of  the  operator  everywhere  in 
the  flow  domain.  By  definition,  H  is  equal  to  one  everywhere  outside  of  a  neighboring  rotor  blade 
row;  at  a  point  within  a  neighboring  rotor  blade  row,  H  is  equal  to  one  for  all  times  that  the  point 
is  immersed  in  fluid  and  equal  to  zero  whenever  the  point  lies  within  a  passing  blade.  The  ratio  of 
time  for  which  the  gate  function  equals  one  to  the  time  for  which  it  equals  zero  is  a  measure  of  the 
geometric  blockage,  \R,  of  the  neighboring  rotor  blade  row. 

The  ensemble-averaged  function,  / ,  can  now  be  decomposed  into  two  components  such 


that 


/-/+/. 


(2.8) 
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where 


f  =  time-averaged  component, 


f  =  deterministic  unsteady  component, 


and,  by  construction, 


/  =  0  . 


Unlike  the  ensemble  averaging  operator,  time  averaging  does  not,  in  general,  commute  with 
differentiation.  Due  to  the  presence  of  the  gate  function,  there  arise  additional  terms  in  the 
interchange  of  averaging  and  differentiation.  The  consequences  of  this  are  discussed  subsequently. 

Time  averaging  the  ensemble-averaged  Navier-Stokes  equation  proceeds  by  decomposing 
the  dependent  variable  in  the  RANS  equations  according  to  Equation  2.8  and  operating  on  the 
resulting  equations  using  Equation  2.7.  In  a  manner  completely  analogous  to  ensemble  averaging, 
the  time  averaging  procedure  gives  rise  to  terms  akin  to  Reynolds  stresses.  These  apparent  or 
"mixing"  stresses  are  a  direct  result  of  the  deterministic  unsteadiness  (that  is,  unsteadiness  that 
correlates  with  shaft  speed)  within  the  flow  field  due  to  the  presence  of  both  rotating  and  stationary 
blade  rows.  Additionally,  however,  the  interchange  of  time  averaging  and  differentiation  gives  rise 
to  forcing  terms  proportional  to  the  average  pressure  loading  and  shear  stress  on  the  surfaces  of  the 
rotating  blades.  Time  averaging  thus  gives  rise  to  a  set  of  body  forces  as  well  as  an  additional  set 
of  mixing  stresses  (or  temporal  correlations)  that  must  be  modelled  along  with  the  Reynolds 
stresses. 

It  is  important  to  recall  the  working  definition  of  the  average-passage  flow  field  for  a  given 
blade  row  embedded  in  a  multiple-blade-row  environment.  The  average-passage  flow  field  for  a 
given  blade  row  is  constructed  to  be  steady  and  spatially  periodic  over  the  pitch  of  the  blade  row. 
For  a  single-stage  machine  (or  a  multi-stage  machine  whose  respective  stator  and  rotor  blade  rows 
all  have  integral-multiple  blade  counts),  subject  to  an  axisymmetric  inflow  and  outflow,  the  time- 
averaged  flow  field  is  the  average-passage  flow  field.  That  is,  ensemble  and  time  averaging  are 
sufficient  to  render  the  flow  field  steady  and  spatially  periodic  over  the  pitch  of  each  of  the  blade 
rows.  Thus,  to  calculate  the  average-passage  flow  for  a  case  such  as  this  requires  only  the 
modelling  of  the  time-averaged  Reynolds  stresses  (from  ensemble  averaging)  and  body  forces  and 
temporal  correlations  (from  time  averaging). 
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2.2.3  Passage-to-Passage  Averaging 

In  more  general  multiple-blade-row  configurations,  the  number  of  blades  in  a  given  blade 
row  is  typically  chosen  such  that  it  is  not  an  integral  multiple  of  the  number  of  blades  in  any  other 
blade  row.  As  such,  in  general,  the  time-average  flow  relative  to  some  blade  row  will  not  be 
spatially  periodic  over  the  pitch  of  that  blade  row.  To  force  this  spatially  aperiodic  flow  field  to  be 
periodic  requires  the  application  of  a  third  averaging  operator.  This  is  the  so-called  passage-to- 
passage  averaging  operator.  Using  the  methods  of  Fourier  analysis,  this  operator  averages  out  the 
passage-to-passage  variations  of  the  time-averaged  flow  field  for  a  given  blade  row  while 
accounting  for  the  global  effect  of  these  variations  through  an  additional  set  of  body  forces  and 
momentum  mixing  correlations.  The  form  of  this  operator  acting  on  the  time-averaged  variable, 


/ ,  is  given  by 
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and 

m  =  number  of  blades  in  the  referenced  blade  row 

Xs  =  geometric  blockage  of  the  neighboring  stationary  blade  row 

G(6)  =  "gate"  function. 

Once  again,  the  purpose  of  the  gate  function  is  to  permit  the  application  of  the  operator 
everywhere  in  the  flow  field.  Analogous  to  the  definition  of  the  gate  function,  H,  for  the  time 
averaging  operator,  G  is  equal  to  one  everywhere  outside  of  neighboring  stator  blade  rows;  within 
neighboring  stator  blade  rows,  G  is  equal  to  one  outside  of  the  blades  and  equal  to  zero  within  the 
blades.  The  ratio  of  the  angular  position  for  which  this  gate  function  equals  one  to  the  angular 
position  for  which  it  equals  zero  is  a  measure  of  the  geometric  blockage,  Xs,  of  the  neighboring 


stator  blade  row.  The  time-averaged  function,  /,  can  now  be  decomposed  into  two  components 
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such  that 


/=/+/, 


(2.10) 


where 


/  =  passage-to-passage-averaged  component 


f  =  spatially  aperiodic  component. 


By  construction,  the  passage-to-passage-averaged  component,  /  ,  is  spatially  periodic  over 


the  pitch  of  the  referenced  blade  row.  Once  again,  unlike  the  ensemble  averaging  operator,  the 
passage-to-passage  averaging  operator  does  not,  in  general,  commute  with  differentiation.  As  with 
the  time  averaging  operator,  the  presence  of  the  gate  function,  G,  in  the  passage-to-passage 
operator  gives  rise  to  additional  terms  from  the  interchange  of  averaging  and  differentiation.  Note 
also  that  because  the  passage-to-passage  averaging  operator  is  constructed  with  its  periodicity  equal 
to  that  of  a  referenced  blade  row,  there  will  be  a  unique  average-passage  flow  field  for  each  blade 
row  in  a  multiple-blade-row  turbomachine. 

Finally,  to  reduce  the  average-passage  equations  from  the  time-averaged  equations  of 
motion  requires  the  application  of  the  decomposition,  Equation  2.10,  and  the  passage-to-passage 
averaging  operator,  Equation  2.9,  to  the  time-averaged  equations.  As  with  time  averaging,  this 
procedure  yields  an  additional  set  of  Reynolds  stress-like  terms.  These  terms  represent  mixing 
correlations  arising  from  passage-to-passage  flow  variations  in  the  time-averaged  flow  field. 
Additionally,  the  interchange  of  passage-to-passage  averaging  and  differentiation  gives  rise  to 
additional  forcing  terms  related  to  the  time-averaged  pressure  and  shear  stress  loadings  on  the  blade 
surfaces  of  neighboring  stator  blade  rows.  In  short,  passage-to-passage  averaging  yields  yet 
another  set  of  body  forces  and  mixing  stresses  (spatial  correlations)  to  be  modelled  in  addition  to 
those  resulting  from  ensemble  and  time  averaging. 
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2.2.4  Modifications  to  the  Navier-Stokes  Equations 

The  average-passage  equations  are  essentially  the  Navier-Stokes  equations  with  some 
additional  terms  and  with  the  dependent  variables  understood  to  be  ensemble,  time,  and  passage-to- 
passage  averaged.  There  are  additional  mixing  stress  terms  that  are  due  to  nonlinear  interactions  of 
turbulent  fluctuations,  deterministic  time-dependent  fluctuations,  and  fluctuations  due  to  spatial 
aperiodicities.  In  essence,  the  shear  stress  terms  in  Equations  2.2  through  2.4  can  be  thought  of  as 
being  augmented  in  the  following  way: 


Tij 


RV 


where  i  and  j  take  on  all  combinations  of  z,  r,  and  8,  and 


R 


=  UU-  +  u.  u. 
v  n  *  j 


u  u 

l  J 


(2.11) 


The  first  term  in  Equation  2.11  represents  the  passage-to-passage,  time-averaged  Reynolds  stress; 
the  second  term  represents  the  passage-to-passage-averaged  mixing  stress  due  to  deterministic 
unsteadiness;  and  the  last  term  represents  a  mixing  stress  due  to  the  steady  aperiodic  velocity  field. 
There  are  also  forcing  functions  that  account  for  the  average  effects  of  pressure  and  shear  stress 
loadings  from  neighboring  rotating  and  stationary  blade  rows.  For  the  sake  of  brevity,  the  full 
equations  are  not  listed  here  (see  Adamczyk  [1984]).  The  modelling  of  the  tensor  (2.11)  as  well  as 
the  body  force  terms  from  time  and  passage-to-passage  averaging  constitutes  the  closure  problem 
for  the  average-passage  system  of  equations. 

In  this  work,  only  the  closure  problem  associated  with  the  ensemble  and  time  averaging  is 
addressed.  As  such,  it  is  more  convenient  to  handle  only  the  Reynolds  stresses  as  additional, 
apparent  stresses  and  to  "lump"  together  all  of  the  mixing  stresses  and  body  forces  from  time 
averaging  into  a  single  source  term.  A  model  for  computing  this  source  is  presented  in 
Section  2.3.  With  this  in  mind,  the  vector  form  of  the  average-passage  equations  for  a  given  blade 
row  rotating  at  a  shaft  speed,  Q  ,  can  be  written  as 

JL(X0  +  il-(XF)  +  i  (XrG)  +  =\K  -\S  ,  (2.12) 

dtx  dz  r  dr  r  36 

where 

F  =  F,-Fv  ,  G  =  G. .  -  Gv,  H  =  Ht-  Hv,  K  =  K.  -  Kv  , 
and  the  subscripts  and  "v"  denote  the  inviscid  and  viscous  portions  of  the  flux  vectors. 
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Equation  2.12  includes  the  following  vector  definitions: 
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where 

we  m  u„  -  rQ 

X  =  geometric  blockage  of  neighboring  blade  rows 

5  =  source  term  containing  the  body  forces  and  correlations  unique  to  the  average  passage 
system  of  equations. 

Several  things  need  to  be  clarified  about  the  above  equations.  First,  note  that  in  each  term  of 
Equation  2.12,  there  appears  a  scale  factor  multiplying  the  dependent  variables.  This  blockage 
factor  is  a  purely  geometrically-derived  scalar.  Relative  to  a  given  blade  row,  it  represents  the 
physical  blockage  due  to  the  presence  of  each  neighboring  blade  row  that  has  been,  in  effect, 
"smeared"  by  the  application  of  the  averaging  operators.  It  is  not  an  additional  dependent  variable; 
it  may  be  calculated  a  priori  for  each  blade  row  in  a  multiple-blade-row  turbomachine. 
Additionally,  the  shear  stresses  in  the  above  viscous  vectors  include  the  additional  Reynolds 
stresses,  that  is,  the  first  term  in  Equation  2.11.  It  is  important  to  note  also  that  in  Equation  2.12, 
all  of  the  dependent  variables  are  interpreted  as  average-passage  variables  (that  is,  the  triple 
overbar  notation  is  assumed);  for  example, 


Tjf,  =*  "rti  ,  et  cetera. 
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In  this  work,  the  integral  form  of  Equation  2.12  is  taken  as  the  basic  system  of  governing 
equations  for  the  average-passage  flow  in  a  multiple-blade-row  environment.  That  is, 

Equation  2.12  is  integrated  over  a  volume  and,  following  an  application  of  Gauss’  theorem,  the 
result  is 


XSdV  .  (2.13) 


Equation  2.13  forms  the  basis  for  the  numerical  treatment  of  the  average-passage  equations. 

In  summary,  carrying  out  the  averaging  procedures  described  in  this  chapter  on  the  Navier- 
Stokes  equations  relative  to  a  given  blade  row,  has  resulted  in  the  following.  As  a  result  of 
ensemble  averaging,  the  random  fluctuations  of  the  dependent  variables  associated  with  turbulent 
flow  have  been  replaced  with  an  additional  mixing  (Reynolds)  stress  acting  throughout  the  flow 
domain.  As  a  result  of  time  averaging,  any  blade  rows  that  are  rotating  relative  to  the  reference 
blade  row  have  been,  in  effect,  "smeared";  their  presence  accounted  for  through  a  body  force 
distribution  in  the  region  occupied  by  the  rotating  blade  and  a  mixing  stress  due  to  the  deterministic 
unsteady  velocity  field  of  the  rotating  blade  row  acting  throughout  the  flow  domain.  As  a  result  of 
passage-to-passage  averaging,  all  blade  rows  that  are  stationary  relative  to  the  reference  blade  row 
and  that  do  not  have  integral-multiple  blade  counts  have  also  been  "smeared";  their  presence 
accounted  for  by  an  additional  body  force  distribution  and  mixing  stress  due  to  the  spatially 
aperiodic  flow  associated  with  the  smeared  blade  row.  The  resulting  average-passage  equations 
thus  consist  of  the  Navier-Stokes  operator  (acting  on  the  average-passage  dependent  variables 
relative  to  a  given  blade  row)  plus  a  series  of  body  force  and  mixing  stress  terms.  The  problem  of 
how  to  model  the  body  force  and  mixing  stress  terms  is  the  closure  problem  for  the  average- 
passage  formulation.  A  simple  closure  model  developed  for  inviscid  flow  through  a  single-stage 
machine  is  briefly  described  in  the  following  section.  Though  in  significantly  less  detail,  it  is  taken 
directly  from  Adamczyk,  Mulac,  and  Celestina  [1986], 


2.3  Closure  Modelling  for  the  Time-Averaged  Equations 

The  focus  of  this  work  is  to  numerically  approximate  a  steady  solution  to  Equation  2.13  for 
a  single-stage  turbomachine.  A  single-stage  machine  consists  of  one  stationary  and  one  rotating 
blade  row.  As  mentioned  previously,  the  time-averaged  equations  in  such  a  case  are  steady  and 
spatially  periodic  and,  as  such,  are  the  average-passage  equations.  The  closure  problem  is  then  to 
model  the  Reynolds  stresses  and  the  body  forces  and  temporal  correlations  associated  with  time¬ 
averaging.  As  such,  the  mixing  stress  takes  the  following  form: 


(2.14) 


=  u<  uj  +  ui  «,  > 


where  the  first  term  is  the  time-averaged  Reynolds  stress  and  the  second  term  is  the  mixing  stress 
due  to  deterministic  unsteadiness. 
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The  Reynolds  stresses  are  grouped  with  the  laminar  stresses  so  that  they  may  be  modelled 
using  an  effective  or  "eddy"  viscosity  construction.  In  this  work,  the  eddy  viscosity  model  is  the 
familiar  algebraic  model  of  Baldwin  and  Lomax  [1978].  It  is  briefly  described  in  the  next  chapter. 

Recall  that  in  Equation  2.13,  all  of  the  average-passage  mixing  stresses  and  body  force 
terms  were  grouped  into  the  source  term,  S.  Adamczyk,  Mulac,  and  Celestina  [1986]  devised  a 
method  for  computing  S  for  an  inviscid  single-stage  machine.  The  method  has  the  very  attractive 
feature  that  the  source  term  for,  say,  blade  row  1  can  be  computed  from  the  average-passage 
solution  for  blade  row  2  and  vice  versa. 

For  the  body  force  calculation,  it  may  easily  be  shown  that  if  it  is  assumed  that  the 
ensemble-  and  time-averaged  pressure  loadings  of  the  blades  are  approximately  equal,  then, 
knowing  the  time-averaged  solution  for  blade  row  1 ,  one  can  calculate  its  body  force  representation 
in  the  frame  of  reference  of  blade  row  2.  This  type  of  body  force  representation  is  axisymmetric. 
Obviously,  this  also  applies  to  the  body  forces  due  to  blade  row  2  acting  in  the  frame  of  reference 
of  blade  row  1 . 

To  model  the  remaining  mixing  stress  terms  (that  is,  the  second  term  in  Equation  2.14)  in 
the  frame  of  reference  of  blade  row  1,  for  example,  the  following  methodology  is  used. 

Decompose  the  i,h  velocity  component  in  the  frame  of  reference  of  blade  row  1  into  two 
components  (for  convenience  the  overbar-hat  notation  is  understood);  that  is, 


where 

m,  =  ensemble-averaged,  deterministic  unsteady  velocity  component  for  blade 
row  1 

m,(2)  =  velocity  component  that  is  steady  with  respect  to  blade  row  2 
u.'  =  velocity  component  that  is  unsteady  with  respect  to  both  blade  rows. 


The  second  term  in  Equation  2.14  then  becomes 


uh 


uPTUjIS  +  u.'u 


m 


«7«r 


As  argued  by  Adamczyk,  Mulac,  and  Celestina  [1986],  the  correlations  involving  u.'  can  be 

expected,  in  general,  to  be  much  smaller  than  those  associated  with  u-2) .  With  this  in  mind,  the 
last  three  terms  in  Equation  2.14  are  neglected  and  the  following  approximation  is  made: 

R9  =  HR?  .  (2-15) 

In  effect,  the  above  approximation  asserts  that  the  dominant  unsteady  correlation  for  blade  row  1  is 
that  due  to  the  steady  hydrodynamic  loading  of  blade  row  2  and  all  correlations  associated  with  the 
blade  row  interaction  velocity  field  are  assumed  to  be  comparatively  small  and  are  neglected. 

Along  with  the  assumptions  regarding  the  body  force  calculation  for  blade  row  1 ,  the 
temporal  correlation,  Equation  2.15,  can  be  shown  to  be  axisymmetric  and  can  be  evaluated  from 
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the  steady  average-passage  flow  field  of  blade  row  2.  Applying  the  same  assumptions  to  the  body 
force  and  temporal  correlations  from  blade  row  1  in  the  frame  of  reference  of  blade  row  2  will 
yield  completely  analogous  results.  As  a  result,  a  rather  simple  equation  can  be  derived  for  the 
total  average-passage  source  for  a  single-stage  turbomachine.  For  the  i,h  blade  row,  it  is  given  by 

| S(Q»)AdV  =  a[  L(0»)  -  jWVv}  -  {  LM(A^)  -  \K{A&»)AdV  }  ,  <2'16) 


where 

L  =  surface  integral  on  the  left-hand  side  of  Equation  2.13 
A  =  circumferential  averaging  operator 
LM  =  axisymmetric  form  of  L 

gw  =  vector  of  average-passage  dependent  variables  for  the  Ith  blade  row. 


In  words,  the  average-passage  source  term  for  a  given  blade  row  is  simply  the  difference 
between  the  circumferential  average  of  the  three-dimensional  flux  operator  acting  on  the  three- 
dimensional  flow  variables  and  the  axisymmetric  flux  operator  acting  on  the  circumferentially- 
averaged,  three-dimensional  flow  variables.  Thus,  the  governing  equations  for  blade  row  1,  for 
example,  will  be  given  by 

|  JL(\&")dV  +  ^[x/*1^  +  XG(1,<£4r  +  XH(')dAe 

=  jXtf'W-  Jxs^v,  (2,17) 

where  the  superscript  refers  to  the  blade  row  and  the  last  term  on  the  right  hand  side  is  given  by 
Equation  2.16  with  i  =  2  (an  analogous  equation  exists  for  blade  row  2). 

Another  result  of  the  assumed  form  of  the  average-passage  source  term  is  the  property  that, 
at  steady-state,  the  circumferential  average  of  the  three-dimensional,  average-passage  flow  fields 
for  both  of  the  blade  rows  will  be  identical.  In  the  numerical  solution  of  the  average-passage 
equations,  this  property  of  the  governing  equations  will  be  useful  as  a  global  convergence  criterion. 


2.4  Pseudo-Compressibility  and  Preconditioning 

The  differential  form  of  the  governing  equations  (Equations  2.12)  does  not  represent  a 
complete  system  of  time  evolution  equations-note  the  zero  element  in  Q.  Because  the  interest  here 
is  in  a  steady-state  solution  to  Equation  2.12,  we  are  free  to  alter  the  time-dependent  portion  of  the 
governing  equation  to  hasten  the  achievement  of  a  steady  solution.  The  basic  idea  is  to  cast 
Equation  2.12  in  a  slightly  modified  form  such  that  the  mathematical  character  of  the  modified 
equation  mimics  that  of  the  compressible  Euler  and  Navier-Stokes  equations.  The  modified 
equations  may  then  be  treated  numerically  using  techniques  developed  for  the  compressible  flow 
equations  of  motion.  This  reformulation  is  referred  to  as  the  "pseudo-compressibility"  approach 
and  was  originally  proposed  by  Chorin  [1967]  and  expanded  upon  by  Turkel  [1986],  among  others. 
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The  pseudo-compressible  form  of  Equation  2.12  is  obtained  by  simply  replacing  Q  in  that 
equation  by  U  where 


P 


U  = 


rue 


In  this  work,  the  time  derivative  in  the  governing  equations  is  multiplied  by  a  "preconditioning" 
matrix,  P,  where 

0  0  0 

1  0  0 

0  1  0 

(cc+l)rug  +  yr2Q  ft  ft  1 

01 

and  a  ,  p  ,  and  y  are  constant  parameters  that  are  used  to  "tune"  the  transient  behavior  of  the 
system  to  accelerate  convergence  to  a  steady-state  solution. 

The  pseudo-compressible  form  of  Equation  2.12  is  thus  given  by  the  following: 

P  JL  (KU)  +  i.  (XF)  +  ii.  (XrG)  +  1A  (Xfl)  =  \K  -  \S  .  (2.18) 

at{  dz  r  dr  r  dd 

Briefly,  some  elaboration  is  needed  regarding  the  choice  of  the  preconditioning  parameters.  It  is 
apparent  that  for  a  =  -1  and  y  =  0  (or  D  =  0),  the  matrix  P  becomes  the  original  "artificial 
compressibility"  preconditioner  developed  by  Chorin  [1967],  with  P  acting  as  a  pseudo-acoustic 
speed.  Also,  for  y  =  0  (or  Q  =  0),  P  becomes  the  preconditioning  matrix  developed  by 
Turkel  [1986],  with  a  and  p  acting  as  parameters  to  "tune"  the  wave  speeds  of  the  system  to 
maximize  convergence  to  steady  state.  Turkel  [1986]  presented  some  guidelines  for  choosing  the 
optimal  values  of  a  and  P  .  In  his  work,  Turkel  [1986]  chose  a  to  be  a  constant  and  P  to  be 
scaled  locally  with  the  velocity  magnitude.  In  this  work,  we  chose  both  a  and  p  to  be  constants. 
Additionally,  the  inclusion  of  the  y  term  in  the  matrix  is  intended  to  render  the  transient  behavior 
of  the  system  independent  of  the  shaft  rotation  rate  fi  .  It  turns  out  that  if  y  is  taken  to  be  equal 
to  -1,  the  eigenvalues  of  the  inviscid  flux  Jacobians  (that  is,  the  signal  propagation  speeds  of  the 
system)  can  be  shown  to  be  algebraically  independent  of  Q  .  Based  on  results  obtained  in  two 
dimensions  for  moving  cascades,  we  decided  to  "hard  wire"  this  value  of  y  for  all  of  the  three- 


P1 

(a+l)Mz 

~W~ 

(a+l)Mr 

~F~ 
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dimensional  calculations.  In  the  cascade  calculations,  we  found  that  for  a  given  choice  of  a  and 
P  ,  the  convergence  to  steady-state  could  be  made  for  any  value  of  the  blade  speed  if  y  was  equal 
to  -1.  For  other  choices  of  y  ,  a  and/or  p  would  sometimes  require  adjustment  to  obtain 
convergence  for  different  blade  speeds.  In  summary,  for  the  numerical  simulation  presented  in 
Chapter  4,  we  used  the  preconditioning  parameters  of  a  =  0.0,  P  =  2.5,  and  y  =  -1.0. 

To  yield  the  equations  that  are  to  be  solved  numerically,  Equation  2.18  is  integrated  over  a 
volume.  The  result  can  be  written  as 

|  pl-{\U)dV  +  <f[\FdL4z  +  \GdAr  +  \HdAe]  =  J  \KdV  -  j  \SdV  .  (2.19) 

It  can  easily  be  shown  that  the  introduction  of  pseudo  compressibility  (that  is,  substituting  the 
vector  U  and  the  matrix  P  into  the  governing  system  of  Equation  2.12)  changes  the  inviscid  form 
of  the  equations  (Re  -*  oo)  to  a  system  of  hyperbolic  equations  analogous  to  the  equations 
governing  inviscid,  compressible  flow.  This  is  useful  in  that  it  permits  the  use  of  the  large  variety 
of  numerical  techniques  developed  for  the  steady  solution  of  the  compressible  equations  of  motion. 
Here,  an  explicit  finite-volume  technique  is  utilized  for  the  numerical  approximation  of  the  steady- 
state  solution  of  Equation  2.19.  This  is  discussed  in  the  following  chapter. 
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3.  Numerical  Discretization 


In  this  chapter,  the  numerical  technique  utilized  to  solve  the  steady-state  average-passage 
equations  for  a  single-stage  turbomachine  is  outlined.  First,  the  spatial  discretization  is  described. 

It  consists  of  a  finite-volume  technique  along  with  a  model  for  an  artificial  dissipation  operator 
necessary  for  numerical  stability.  This  is  followed  with  a  description  of  the  explicit  multistage  time 
integration  procedure  utilized  along  with  several  convergence  acceleration  techniques-local  time 
stepping,  implicit  residual  averaging,  and  multigrid.  Next,  the  boundary  conditions  required  by  the 
average-passage  simulation  are  briefly  described.  This  chapter  ends  with  a  description  of  the 
numerical  treatment  of  the  average-passage  closure  models. 


3.1  Cell-Centered  Finite-Volume  Discretization  in  Space 

A  generic  flow  domain  is  schematically  depicted  in  Figure  3.1  where  the  relevant 
coordinate  systems  are  identified.  This  domain  is  discretized  into  a  series  of  contiguous  hexahedral 
cells.  The  numerical  approximation  of  the  governing  system  of  equations,  Equation  2.19,  begins 
with  the  standard  cell-centered  finite- volume  discretization  for  hexahedra.  That  is,  for  a  hexahedral 
cell,  Equation  2.19  is  numerically  integrated  in  space  to  yield 

P*(KV)'V«m  *  £[<XR4,)„  *  (XG A).  *  (XftV J„ 

at\  «->  (3.1) 

=  <x*>*w*  -  , 

where  the  subscript  ”ijk”  denotes  the  (i,j,k)th  control  volume  or  cell,  m  denotes  the  m'h  face  of  the 
cell,  V ol  denotes  the  cell  volume,  and  Av  Ar,  and  Ag  are  components  of  the  outward  directed  area 
of  a  cell  face. 

In  Equation  3.1,  the  dependent  variables  are  defined  at  the  center  of  a  control  volume  and 
are  taken  to  be  an  average  value  for  the  cell.  Figure  3.2  depicts  the  (ij,kf  cell  with  its  defining 
grid  points.  The  cell  volume  is  computed  using  a  tetrahedron  decomposition  procedure.  The 
directed  areas  in  Equation  3.1  are  computed  using  the  cross  product  of  the  diagonals  of  each  cell 
face.  Figure  3.3  shows  the  directed  areas  of  a  typical  cell. 

If  all  of  the  spatial  terms  in  Equation  3.1  are  collectively  denoted  by  E,  then  Equation  3.1 
may  be  written  more  simply  as 


*  f3W|.  -  0  •  (3-2) 

The  portions  of  E  representing  facial  fluxes  are  approximated  using  simple  arithmetic  averaging  of 
neighboring  cell-centered  quantities.  As  such,  it  may  be  easily  shown  that  for  uniformly  spaced 
grids  in  cylindrical  coordinates,  Equation  3.2,  in  its  present  form,  amounts  to  a  central-difference 
approximation  to  Equation  2.18.  It  is  well  known  that  central -difference  approximations  of  the 
Euler  equations  or  high  Reynolds  number  Navier-Stokes  equations  requires  the  addition  of  some 
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artificial  damping  terms  for  stability.  With  the  inclusion  of  these  terms,  Equation  3.2  becomes 


(3.3) 


where  D  is  the  artificial  dissipation  operator. 

3.2  Artificial  Dissipation 

The  artificial  dissipation  operator  used  in  this  work  is  modelled  after  the  operator  originally 
proposed  by  Jameson,  Schmidt,  and  Turkel  [1981]  for  transonic  flow,  with  modifications  for 
highly-stretched  grids  suggested  by  Martinelli  [1987].  In  the  present  implementation,  the 
assumption  of  pseudo-compressibility  results  in  some  slight  modifications  and  simplifications. 

The  purpose  of  adding  artificial  dissipation  terms  is  to  provide  damping  of  error  waves  in 
the  solution  domain.  An  effective  operator  can  be  constructed  by  looking  at  the  natural  dissipation 
inherent  in  upwind  differencing.  In  other  words,  a  form  for  the  dissipation  operator  can  be 
suggested  by  determining  what  form  D  should  take  so  that  when  it  is  combined  with  the  central- 
differenced,  preconditioned  convection  terms  (as  in  Equation  3.3),  the  result  will  approximate  an 
upwind  differencing  of  just  the  preconditioned  convection  terms.  The  actual  form  of  the 
dissipation  will  depend  on  the  type  of  upwind  differencing  considered;  first-order  upwinding  will 
result  in  a  Laplacian  form  for  D,  while  second-order  upwinding  will  yield  a  D  of  biharmonic  form. 
In  both  cases,  the  proper  form  for  the  dissipation  will  scale  with  the  signal  propagation  speeds  (or 
wave  speeds)  of  the  preconditioned  convection  terms.  The  signal  propagation  speeds  are 
determined  by  the  eigenvalues  of  the  flux  Jacobian  matrices  formed  by  the  inviscid  flux  vectors  Ft, 


G„  and  Ht. 


Jameson,  Schmidt,  and  Turkel  [1981]  originally  devised  an  artificial  dissipation  operator 
that  was  a  combination  of  second-  and  fourth-difference  operators  with  the  relative  amounts  of  each 
determined  by  the  flow  solution.  Additionally,  both  terms  in  his  operator  were  scaled  by  the 
maximum  wave  speeds  or,  more  precisely,  by  the  spectral  radii  of  the  inviscid  flux  Jacobians. 
Although  their  original  implementation  was  for  transonic  flow,  the  same  methodology  for  the 
construction  of  the  dissipation  operator  is  followed  in  this  work. 

In  a  three-dimensional  application  (such  as  in  this  work),  the  dissipation  operator  is 
constructed  separately  in  each  of  the  coordinate  directions  using  the  above  methodology.  The 
general  form  of  the  operator  can  thus  be  written  as 


(3.4) 


where  ri ,  C ,  and  i  denote  the  i,  j,  and  k  computational  directions,  respectively.  Note,  in  the 
following,  for  illustrative  purposes,  only  the  ti -direction  is  considered  in  detail;  the  other 
directions  follow  analogously. 
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Following  the  original  formulation,  a  fully  conservative  dissipation  operator  is  constructed 
containing  both  second  and  fourth  differences  with,  for  example, 


W  = 


(3.5a) 


where 


=  <«?. 


A  - 

V 


.(4) 


m 


1 


(3.5b) 


and  and  A^  are  undivided  first  and  third  difference  operators  in  the  q  -direction  (note  that 

the  subscripts  j  and  k  are  implied  in  all  terms).  For  the  coefficients  e(2)  and  e(4),  eigenvalue 
scaling  is  used.  Specifically,  these  coefficients  are  defined  to  be  proportional  to  a  blending  of  the 
spectral  radii  of  the  preconditioned  inviscid  flux  Jacobians  in  Equation  3.3;  the  reason  for  the 
blending  will  be  described  shortly.  The  second  and  fourth-difference  dissipation  scale  factors  are 
given  by 


e 


(2) 


—  k(2)<£  ,  min(  v  .  i 

2  * 


,1.0) 


(3.6a) 


and 


where 


e<4  ,  =  max(i- 
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,0.0) 


(3.6b) 


v,.l  =  raax(v.,vM) 
1 


-  2_Pi  +_  Pi o 

+  2  Pi  +  Pm 


(3.7) 
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In  the  above,  ,  i|ff ,  and  i|i5  are  the  maximum  wave  speeds  in  the  respective  computational 
coordinate  directions  and  A  is  the  geometric  blockage  factor.  For  example,  in  the  -direction, 


YV 


(3.9) 


where 


<?,  =  (!-  ccKV-A")  , 

5  =  A”  -A”  , 

Av  =  Alet  +  A'rer  +  AJet  , 


and 

V  =  velocity  vector 

A1  =  the  directed  area  of  an  r|  =  constant  cell  face 

ez,  er,  and  eQ  =  unit  vectors  in  thez,  r,  and  0  directions,  respectively. 

Also  in  the  above,  (J)^  is  the  blending  function  for  the  eigenvalue  scaling  in  each  of  the 
coordinate  directions.  This  blending  has  been  found  to  be  useful  for  flow  solutions  using  highly- 
stretched  computational  grids  by  Martinelli  [1987].  The  purpose  of  this  function  is  to  combine  the 
artificial  dissipation  scale  factors  (that  is,  the  spectral  radii  of  the  inviscid  flux  Jacobians)  in  each  of 
the  coordinate  directions  so  that  no  one  coordinate  direction  is  artificially  damped  significantly 
more  than  any  other.  Note  that  the  exponent  6  in  the  definition  of  the  blending  function 
determines  the  level  of  blending;  for  example, 

6=0  scales  the  dissipation  independently  in  each  direction  (anisotropically) 


and 


6=1  scales  the  dissipation  identically  in  each  direction  (isotropically). 

For  typical  high  Reynolds  number  grids,  6  «  0.50  has  been  found  to  be  adequate. 

There  are  two  other  constants  that  have  to  be  set  by  the  user  in  Equation  3.5b;  these  are 
k(2)  and  k(4).  The  factors  of  1/2  and  1/4  in  Equations  3.6a  and  3.6b  are  included  so  as  to  confine 
the  typical  values  of  k®  and  k(4)  to  the  range  of  0.0  to  1.0. 

In  Equation  (3.6a),  v  can  be  interpreted  as  a  pressure  gradient  "switch"  that  locally  alters 
the  artificial  dissipation  depending  on  the  gradients  present  in  the  flow  field.  This  function 
essentially  measures  the  rate  of  change  of  the  pressure  gradient  in  a  particular  coordinate  direction. 
The  way  in  which  this  function  is  incorporated  in  the  definitions  of  the  artificial  dissipation  scale 
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factors  results  in  e(4)  tending  to  zero  in  regions  of  rapid  change  in  pressure  gradient;  while  in  the 

same  regions,  e(2)  approaches  its  maximum  value.  Conversely,  in  regions  where  the  pressure 
gradient  is  only  slowly  varying,  the  second-difference  artificial  dissipation  will  approach  zero  and 
there  will  remain  only  a  background  level  of  fourth-difference  artificial  dissipation.  In  the  original 
(compressible  flow)  implementation  of  this  type  of  blended  artificial  dissipation,  the  intention  was 
the  sharp  capture  of  shock  waves.  Near  shock  waves  it  was  found  that  the  fourth-difference 
artificial  dissipation  would  cause  overshoots  that  could  negatively  effect  convergence.  If,  however, 
the  fourth-difference  portion  of  the  operator  was  turned  off  in  the  vicinity  of  a  shock  and  the 
second-difference  portion  increased,  shocks  could  be  captured  without  overshoots.  Of  course,  this 
rendered  the  solution  only  first-order  accurate  near  shocks,  but  it  would  remain  second-order 
accurate  in  smooth  regions  of  the  flow  field.  Although  this  work  deals  exclusively  with 
incompressible,  shock-free  flow,  this  type  of  formulation  is  preserved  for  handling  the  possibility 
of  nearly  discontinuous,  shock-like  structures  that  may  appear  during  the  transient.  In  most 
situations,  however,  this  blended  dissipation  is  not  necessary  for  convergence.  For  this  reason,  the 
following  special  cases  are  included  in  the  formulation: 

If  k(4)  =  0,  then 


and 


_(2) 


—  K®  <j>  , 

2  ” 


0.0 


while  if  k(2)  =  0,  then 


e(2) ,  =  0.0 

t+- 

2 


and 


In  other  words,  by  choosing  either  of  the  constants  to  be  zero,  the  dissipation  operator  can  be  made 
independent  of  the  switching  function  and,  as  such,  represent  a  pure  second-difference  or  a  pure 

fourth-difference  dissipation.  Therefore,  in  these  special  cases,  the  coefficients  k(2)  and  k(4) 
approximately  represent  the  fraction  of  the  equivalent  dissipation  that  would  be  present  from  first- 
order  or  second-order  upwind  discretizations  of  the  inviscid  convection  terms  in  the  governing 

equations.  For  example,  choosing  k(2)  =  0.10  and  k(4)  =  0.0  implies  a  level  of  numerical 
dissipation  of  approximately  ten  percent  of  the  level  that  would  be  present  if  first-order  upwinding 
were  used  for  the  convection  terms. 

The  boundary  treatment  for  the  artificial  dissipation  follows  that  of  Rizzi  and 
Eriksson  [1985].  Basically,  the  boundary  and  near-boundary  cells  are  handled  in  a  way  that 
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ensures  that  both  the  second-  and  fourth-difference  operators  remain  dissipative  at  the  extremes  of 
the  domain. 


3.3  Time  Integration 

To  facilitate  the  description  of  the  time  stepping  methodology,  the  residual  is  defined  as 

K(tO„  -  PlE(U)„  -  D(U)„  ■  <3'10> 

It  is  important  to  note  that  this  residual  is  composed  of  a  number  of  distinct  elements— recall  that 
the  vector  E  contains  all  of  the  discretized,  physical  spatial  terms  (that  is,  inviscid  and  viscous 
fluxes,  as  well  as  the  radial  momentum  and  average-passage  source  terms)  and  the  vector  D  is  the 
artificial  dissipation  operator  described  in  Section  3.2.  This  has  important  consequences  in  the 
implementation  of  the  time  stepping  scheme.  Using  the  residual  definition,  Equation  3.10,  the 
discrete  governing  equations  may  be  written  as 

A.(mf,Voltt  *  R(U)fl  -  0  .  (3.11) 

Equation  3.11  is  then  integrated  in  time  to  a  steady-state  solution  using  the  explicit  "hybrid 
multistage"  scheme  of  Jameson,  Schmidt,  and  Turkel  [1981].  In  its  simplest  manifestation,  this 
four-stage  Runge-Kutta-like  scheme  takes  on  the  following  form  for  the  ( i,j,k)'h  cell  at  the  n'h  time 
step: 


Um  =  Un  -  a,  -^L  R(Un) 

1  Wol  ' 

U{2)  =  Un  -  a 2-^LR(W") 

2  Wol 

UO)  =  un  -  a,  ^  R(U2)) 

3  Wol 

Un"  =  Un  - 

4  Wol 


(3.12) 


where 

al  =  0.250,  a2  =  0.333,  a3  =  0.500,  a4  =  1.000  , 


and  the  bracketed  numerical  superscripts  denote  provisional  values. 

It  turns  out  that  evaluating  the  entire  residual  at  each  stage  in  Equation  3.12  can  be  very 
expensive  (for  example,  the  artificial  dissipation  operator,  D,  is  roughly  twice  as  expensive  to 
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evaluate  as  just  the  inviscid  portion  of  E).  In  practice,  it  has  been  found  that  it  is  sufficient  to 
evaluate  the  physical  dissipation  and  artificial  dissipation  portions  of  R  at  the  first  stage  only  and 
freeze  their  values  through  subsequent  stages.  This  technique  has  only  a  very  small  detrimental 
effect  on  the  size  of  the  permissible  time  step,  while  saving  significant  computational  time.  This 
selective  evaluation  of  the  different  pieces  of  the  residual  at  each  stage  is  what  distinguishes  this 
type  of  scheme  from  a  true  Runge-Kutta  scheme. 

In  this  work,  a  slightly  different  four-stage  scheme  is  used.  The  coefficients  for  the 
scheme  presently  used  are 

a,  =  0.250,  a2  =  0.500,  a3  =  0.550,  a4  =  1.000  . 

This  scheme  was  chosen  because  of  its  more  effective  high  wave  number  damping  characteristics 
relative  to  the  "standard"  scheme.  Although  the  permissible  time  step  for  this  scheme  is  slightly 
less  than  the  standard  scheme,  its  high  wave  number  damping  gives  it  superior  multigrid 
performance  and  overall  faster  convergence. 


3.4  Convergence  Acceleration 

Several  techniques  are  utilized  to  hasten  the  attainment  of  a  steady-state  solution  to 
Equation  3.3.  The  techniques  used  are:  local  time  stepping,  implicit  residual  averaging,  and 
multigrid. 


3.4.1  Local  Time  Stepping 

Because  the  time  coordinate  in  the  governing  equations  is  purely  artificial,  the  transient 
behavior  of  the  system  may  be  freely  manipulated  to  accelerate  convergence  to  steady  state  (the 
introduction  of  the  preconditioning  matrix,  P,  was  the  first  instance  of  this  sort  of  manipulation). 
The  rate  of  convergence  is  determined  by  how  quickly  errors  can  be  removed  from  the  solution 
domain.  There  are  essentially  three  mechanisms  by  which  errors  can  be  removed  from  the 
domain.  They  can  be  damped  by  numerical  (or  physical)  diffusion;  they  can  be  transported  with 
the  flow  through  the  outflow  boundary;  or  they  can  propagate  "acoustically"  out  of  the  domain  in 
any  direction.  Excessive  numerical  damping  will  have  a  detrimental  effect  on  the  steady-state 
solution.  As  such,  the  amount  of  artificial  dissipation  is  generally  kept  as  low  as  possible.  The 
most  effective  way  to  increase  the  convergence  rate,  therefore,  is  to  maximize  the  computational 
rate  at  which  information  is  propagated  throughout  the  domain,  or,  in  other  words,  to  minimize  the 
number  of  time  steps  necessary  to  propagate  a  signal  over  a  given  distance.  This  is  most  easily 
done  by  advancing  the  solution  for  each  cell  in  the  domain  through  the  maximum  possible  time  step 
for  that  cell. 

To  determine  the  maximum  allowable  time  step,  a  Fourier  analysis  is  carried  out  on  a 
linearized  form  of  the  governing  equations  in  generalized  curvilinear  coordinates.  The  result  of 
this  analysis  is  a  time  step  determined  by  two  conditions  imposed  by  the  inviscid  and  viscous 
portions  of  the  governing  equations.  The  inviscid  limitation  on  the  permissible  time  step  is  the 
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familiar  Courant-Friedrichs-Lewy  ( CFL )  condition  for  the  hyperbolic  Euler  equations.  In  this  case, 
it  is  formed  directionally  in  the  following  way: 


<*!>* 


CFL 


1  1 

+  _  -f  _ 

A  L  A L 


ijk 


(3.13) 


where  the  subscript  "I"  denotes  inviscid  and  CFL  is  a  stable  Courant  number  for  the  multistage 
scheme  in  Equation  3.12.  The  directional  time  steps  are  functions  of  the  maximum  wave  speeds 
(flux  Jacobian  spectral  radii)  in  the  respective  computational  coordinate  directions.  For  example, 
in  the  q  -direction, 


(3.14) 


where  is  as  given  in  Equation  3.9;  however,  here  the  average  value  for  a  cell  is  used.  The  £ 
and  l  directions  follow  analogously. 

The  viscous  limitation  on  the  permissible  time  step  is  a  result  of  the  parabolic  or  diffusive 
portions  of  the  governing  equations.  As  such,  it  has  terms  dependent  on  the  physical  as  well  as 
artificial  dissipation  operators  in  Equation  3.3.  To  derive  the  viscous  stability  bound,  the 
methodology  described  by  Kunz  and  Lakshimarayana  [1992]  was  applied  to  the  preconditioned 
pseudo-compressible  equations  of  motion.  The  resulting  viscous  time  step  limit  can  be  expressed 
in  the  following  way  if  the  assumption  of  grid  orthogonality  is  made: 
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(3.15) 


where  the  subscript  "V"  denotes  viscous  and  VON  is  a  stable  von  Neumann  number  for  the 
multistage  scheme.  The  directional  factors  in  the  denominator  have  the  following  form: 


<rA*  = 
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where 


VT 


1  + 
Re 


and  ve  is  the  eddy  viscosity.  The  values  of  e®  and  e^4)  for  a  cell  are  computed  using  the  average 

values  of  Equations  3.6a  and  3.6b  in  the  tj  -direction  for  that  cell.  The  £  and  £  directions  follow 
analogously. 
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The  value  of  the  allowable  time  step  for  a  particular  cell  is  then  determined  by  simply 
taking  the  more  restrictive  of  the  two  limits  (Equations  3.13  and  3. 15)— that  is, 

Atijk  =  min(A tnAtv)ijk  .  (3.16) 

Stable  values  for  CFL  and  VON  can  be  determined  from  an  examination  of  the  amplification  factor 
of  the  multistage  scheme  (Equation  3.12).  Figure  3.4  shows  contours  of  the  amplification  factor 
for  the  scheme  used  in  this  work.  The  outer  boundary  of  this  region  is  the  unit  amplification 
contour— that  is,  neutral  stability.  The  maximum  allowable  Courant  number  is  equal  to  the 
intercept  of  the  unit  contour  with  the  imaginary  axis,  while  the  maximum  allowable  von  Neumann 
number  is  given  by  the  unit  contour’s  intersection  with  the  negative  real  axis.  A  safe  way  to 
choose  CFL  and  VON  is  to  choose  them  such  that  the  rectangle  that  they  define  lies  completely 
within  the  stability  region.  For  the  multistage  scheme  used  in  this  work,  a  good  choice  is 

CFL  =  2.0  ,  VON  =  1.0  , 
with  the  corresponding  rectangle  shown  in  Figure  3.4. 


3.4.2  Implicit  Residual  Averaging 

Another  widely  used  technique  for  accelerating  the  convergence  of  Runge-Kutta-type 
schemes  is  residual  averaging  (or  smoothing).  The  basic  idea  is  to  increase  the  allowable  Courant 
and  von  Neumann  numbers  (and,  hence,  the  size  of  the  local  time  step)  by  increasing  the  spatial 
support  of  the  discrete  approximation.  One  way  to  do  this  is  to  increase  the  number  of  stages  in 
the  multistage  scheme.  Unfortunately,  this  is  rather  expensive  computationally.  Another  way  to 
increase  spatial  support  is  to  spatially  average  the  residuals  in  some  way.  That  is,  replace  the 
residual,  R,  for  a  given  cell  with  a  new  residual  computed  using  residual  information  from 
neighboring  cells.  Jameson  [1983]  found  it  most  effective  to  carry  out  this  averaging  implicitly.  It 
turns  out  that  increasing  the  time  step  a  given  amount  in  this  manner  is  cheaper  computationally 
than  adding  more  stages  to  the  time  integration.  Also,  implicit  residual  averaging  has  no  imposed 
upper  bound  on  the  size  of  the  increase  in  time  step,  as  does  the  inclusion  of  additional  stages 
(there  are  practical  limits,  however). 

Basically,  residual  averaging  consists  of  the  replacement  of  the  residual,  R,  for  a  cell, 
defined  by  Equation  3.10,  with  an  implicitly-averaged  residual,  R*,  where  R*  is  the  solution  to  the 
following: 


(l-SA„)(l-6fArr)(l-efA{{)/?*  =  R  .  (3.17) 

In  Equation  3.17,  ,  ACf,  and  A?c  are  undivided  second-difference  operators  and  e  ,  ef,  e? 

are  smoothing  factors  in  the  respective  computational  coordinate  directions.  In  this  work,  the 
smoothing  factors  are  taken  to  be  constants  typically  in  the  range  of  0.0  to  1.0  (For  the  results 
shown  in  Chapter  4  they  were  each  equal  to  1.0.).  Equation  3.17  is  solved  by  three  successive 
block  tridiagonal  sweeps  through  the  domain.  In  practice,  it  is  sufficient  to  replace  the  residuals  in 
Equation  3.12  at  the  second  and  last  stages  only.  This  is  sufficient  to  approximately  double  the 
allowable  time  step. 
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3.4.3  Multigrid 


A  multigrid  strategy  based  on  the  technique  outlined  by  Jameson  and  Baker  [1984]  is  also 
implemented  as  a  convergence  acceleration  technique.  Its  implementation  in  a  compressible 
average-passage  formulation  is  outlined  by  Kirtley,  Beach,  and  Adamczyk  [1990],  Here,  its 
application  to  the  preconditioned,  pseudo-compressible,  average-passage  equations  is  briefly 
described. 


The  primary  advantage  of  using  a  multigrid  technique  for  the  steady-state  solution  of  the 
Euler  or  Navier-Stokes  equations  is  the  effectiveness  of  the  multigrid  procedure  in  eliminating  low 
wave  number  errors  from  the  solution  domain.  This  is  especially  attractive  in  the  current 
application  to  turbomachine  simulation  for  the  following  reason:  The  boundary  conditions  used  in 
this  work  require  that  the  exit  hub  static  pressure  be  periodically  adjusted  (or  throttled)  to  draw  the 
correct  mass  flow  through  the  machine.  In  the  case  of  rotors,  it  is  essential  that  the  correct  (or 
design)  mass  flow  be  simulated,  as  otherwise,  the  flow  incidence  angles  and,  hence,  pressure 
loading  of  the  rotor  blades  will  be  incorrect.  Each  time  the  back  pressure  is  adjusted,  a  pressure 
wave  is  sent  upstream  to  readjust  the  flow  accordingly.  These  pressure  waves  are  of  very  long 
wavelength  and,  as  such,  are  slow  to  be  removed  by  the  multistage  algorithm.  The  multigrid 
procedure  is  a  means  by  which  long  wavelength  errors  like  these  can  be  quickly  removed  by 
performing  intermediate  iterations  on  coarsened  grids  where  the  waves  appear  with  higher  wave 
numbers. 

In  this  work,  a  standard  "V-cycle"  multigrid  algorithm  is  utilized.  The  cycle  consists  of  a 
series  of  multistage  iterations  on  progressively  coarser  grids.  Each  coarse  grid  iteration  includes  a 
forcing  function  derived  from  information  obtained  from  finer  grid  solutions.  Once  the  coarsest 
grid  in  the  sequence  is  reached  and  the  corrections  to  the  flow  variables  on  this  grid  are  computed 
by  advancing  the  solution  through  a  multistage  iteration  (or  time  step),  a  series  of  interpolation 
steps  follow  where  corrections  on  each  coarsened  grid  are  interpolated  to  progressively  finer  grids. 
Ultimately,  the  flow  variables  on  the  finest  grid  are  corrected  and  the  cycle  can  begin  again. 

Although  the  multigrid  algorithm  can  theoretically  be  carried  out  over  as  many  grid 
coarsenings  as  the  finest  grid  can  be  evenly  divided  into,  the  current  implementation  of  the  scheme 
is  limited  to  two  grid  coarsenings-that  is,  a  three-level,  V-cycle  multigrid. 


Specifically,  upon  completion  of  a  fine-grid  iteration  (or  time  step),  the  flow  variables  on 
the  fine  (h)  grid  are  restricted  to  a  coarsened  ( 2h )  grid  by  the  following  volume-weighted  average: 


U, 


2 h 


5>„I \VoIt 
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(3.18) 


where  the  summation  is  over  the  eight  fine-grid  cells  that  make  up  one  coarse-grid  cell.  A  forcing 
function  is  computed  so  that  fine-grid  accuracy  is  maintained  for  an  iteration  on  the  coarsened  grid. 
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For  the  2h  grid,  it  is  given  as 


/»  =  -  Ew, , 

8 


(3.19) 


where  RfU)^  is  the  residual  computed  on  the  coarse  grid  using  the  restricted  variables  from  the 
fine  grid.  With  the  above  definitions,  the  m'h  stage  of  the  multistage  scheme  on  the  2h  grid  will 
look  like 
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(3.20) 


One  further  coarsening  is  carried  out  to  a  4h  or  doubly-coarsened  grid,  with  U4h  defined 
analogously  to  Equation  3.18.  The  forcing  function  on  this  grid  is  given  as 

/«  -TO.  -  £{TO»  -/»}  ■  (3.21) 
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The  m*  stage  of  the  multistage  scheme  for  the  4h  grid  is  then  simply  given  by  Equation  3.20  with 
the  2h  subscripts  everywhere  replaced  by  4h.  Upon  completion  of  a  time  step  on  the  4h  grid,  the 
corrections  to  the  4h  flow  variables  are  then  prolongated  (interpolated)  to  the  2h  grid  using  simple 
tri-linear  interpolation.  The  modified  2h  flow  variables  are  then  used  to  calculate  corrections  to  the 
flow  variables  on  the  2h  grid.  These  corrections  are  then  prolongated  to  the  fine  grid.  The 
corrected  variables  on  the  fine  grid  are  then  used  to  initiate  the  next  fine-grid  iteration  and  a  three- 
level  multigrid  cycle  is  complete. 

It  is  important  to  note  that  on  both  of  the  coarsened  grids  {2h  and  4h),  only  the  inviscid 
portions  of  the  residuals  are  computed.  That  is,  RfU)^  and  R(U)^  are  essentially  Euler  operators 
and,  as  such,  no-flux  solid  surface  velocity  boundary  conditions  are  applied  on  these  grids 
regardless  of  whether  the  fine-grid  flow  is  viscous  or  inviscid.  All  other  boundary  conditions  on 
the  coarsened  grids  are  the  same  as  those  used  on  the  fine  grid.  Also,  just  as  in  the  case  of  a 
single-grid  calculation,  the  coarsened-grid  residual  operators  include  an  artificial  dissipation 
operator.  Unlike  the  fine-grid  calculation,  however,  the  operator  on  2h  and  4h  grids  is  a  pure 
second-difference  operator.  This  is  done  to  ensure  effective  high  wave  number  damping  on  the 
coarsened  grids.  Finally,  all  of  the  viscous  effects  and  multiple-blade-row  effects  (average-passage 
source)  are  "felt"  on  the  coarsened  grids  through  the  presence  of  the  fine-grid  residual  in  the 
definitions  of  the  forcing  functions,  /a  and 


3.5  Boundary  Conditions 

In  a  blade  passage,  there  are  four  distinct  types  of  boundaries:  inflow,  outflow,  solid 
surface,  and  periodic.  Figure  3.5  depicts  schematically  the  boundary  condition  types  and  locations 
for  a  typical  blade  passage.  As  this  is  a  cell-centered  finite-volume  discretization,  all  boundary 
conditions  are  satisfied  with  the  use  of  "phantom"  cells  located  outside  of  the  physical  domain. 

The  boundary  conditions  utilized  for  solid  surfaces  differ  depending  on  whether  the  flow  is  viscous 


27 


or  inviscid.  For  all  other  types  of  boundaries,  the  boundary  conditions  are  the  same  for  either  type 
of  flow. 


Consistent  with  the  hyperbolic  character  of  the  inviscid  preconditioned  pseudo-compressible 
equations  of  motion,  at  an  inflow  boundary,  three  conditions  must  be  specified  and  a  fourth 
extrapolated  from  inside  the  flow  domain.  In  this  work,  for  both  inviscid  as  well  as  viscous  flow, 
the  total  pressure  and  two  absolute  flow  angles  are  specified  at  the  inflow  plane;  the  static  pressure 
is  extrapolated  from  the  interior. 

At  an  outflow  boundary,  the  character  of  the  governing  equations  dictates  that  one 
condition  be  specified  and  three  other  conditions  extrapolated  from  the  interior.  Here,  the  static 
pressure  along  the  hub  at  the  outflow  plane  is  specified  and  integrated  outward  using  simplified 
radial  equilibrium— that  is, 


dp  =  Ue_  (3.22) 

dr  r 

Also  at  this  boundary,  all  three  velocity  components  are  extrapolated  from  the  interior  solution. 

Recall  that  the  governing  preconditioned  pseudo-compressible  average-passage  equations 
are,  by  definition,  periodic  over  the  pitch  of  a  blade  passage.  As  such,  all  flow  variables  are  taken 
to  be  periodic  upstream  and  downstream  of  the  blade  surfaces  and  in  any  hub  or  tip  gap  regions 
above  or  below  the  blades. 

On  the  blade,  hub,  and  shroud  (casing)  surfaces,  different  boundary  conditions  are 
employed  depending  on  whether  the  flow  is  inviscid  or  viscous.  For  inviscid  flow,  the  solid 
surface  boundary  condition  for  velocity  is  the  familiar  no-flux  condition,  while  the  pressure  is 
obtained  by  an  application  of  the  normal  momentum  equation  for  the  surface  in  question.  For 
viscous  flows,  where  there  is  adequate  grid  resolution  in  the  boundary  layer,  a  no-slip  condition  is 
satisfied  for  the  velocity  components.  For  higher  Reynolds  number  calculations  where  grid 
resolution  is  a  problem,  there  is  also  a  viscous  slip  velocity  condition,  where  the  slip  velocity  is 
determined  by  a  law-of-the-wall-based  wall  function.  For  both  cases,  the  wall  pressure  is 
determined  by  assuming  a  zero  normal  pressure  gradient  at  the  wall.  Also,  as  mentioned  by 
Kirtley,  Beach,  and  Adamczyk  [1990],  some  minor  additional  terms  are  included  in  the  hub  and 
shroud  boundary  conditions  to  ensure  that  a  common  axisymmetric  solution  is  converged  upon  by 
all  blade  rows. 


3.6  Closure  Modelling 

This  work  deals  with  the  calculation  of  the  average-passage  flow  fields  associated  with  a 
single  stage  turbomachine— that  is,  two  blade  rows.  As  such,  the  closure  issues  are  concerned  with 
only  the  ensemble  averaging  and  the  time  averaging.  Specifically,  these  are  the  first  two  terms  in 
the  mixing  stress  relation  of  Equation  2.11,  in  addition  to  the  blade  force  terms  resulting  from  the 
time  averaging  of  the  momentum  equations. 
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3.6.1  Reynolds  Stress  Closure 

The  time-averaged  Reynolds  stress  (the  first  term  in  Equation  2.11)  is  modelled  using  a 
familiar  eddy  viscosity  construction.  Specifically,  the  algebraic  eddy  viscosity  model  of  Baldwin 
and  Lomax  [1978]  is  used.  This  model  is  a  two-layer  mixing-length-type  model  where  the  relevant 
length  scale  is  based  on  the  distance  to  the  nearest  wall.  In  the  current  application,  no  special 
corner  treatment  is  utilized  In  the  wakes,  a  simple  extrapolation  procedure  is  used  based  on  the 
eddy  viscosity  distribution  at  the  trailing  edge.  This  extrapolation  procedure  assumes  that  the  axial 
grid  lines  approximately  follow  the  downstream  trajectory  of  the  wake.  The  eddy  viscocity 
distribution  at  the  trailing  edge  is  then  exponentially  decayed  along  these  grid  lines,  as  suggested 
by  Adamczyk  [1992],  To  ensure  that  a  common  circumferentially-averaged  flow  field  is  obtained 
for  all  blade  rows  in  a  multiple-blade-row  environment,  the  computed  eddy  viscosity  for  a  blade 
passage  is  modified  as  described  by  Kirtley,  Beach,  and  Adamczyk  [1990], 


3.6.2  Average-Passage  Closure 

The  temporal  correlations  (mixing  stresses)  and  body  forces  resulting  from  the  application 
of  the  time-averaging  operator  to  the  Reynolds-averaged  Navier-Stokes  equations  are  computed 
using  a  discretized  version  of  the  technique  described  in  Section  2.3.  Recall  that  the  methodology 
in  Section  2.3  was  formulated  to  calculate  only  that  portion  of  the  temporal  correlation  for  a  given 
blade  row  due  to  the  velocity  components  that  are  steady  with  respect  to  neighboring  rotating  blade 
rows.  What  follows  is  a  brief  outline  of  how,  in  practice,  the  average-passage  flow  fields  for  a 
two-blade-row  turbomachine  are  determined.  It  follows  the  method  described  by  Adamczyk, 
Mulac,  and  Celestina  [1986] 

The  form  of  the  governing  average-passage  equations  suggested  by  Equation  2.16  implies 
that  the  flow  fields  are  to  be  calculated  simultaneously.  In  practice,  however,  a  two-level  iteration 
procedure  is  used.  In  the  inner  iterations,  the  flow  fields  for  the  individual  blade  passages  are 
updated  using  the  multistage  algorithm  with  the  average-passage  source  term  fixed.  In  the  outer 
iteration,  the  average-passage  source  terms  are  updated  according  to  the  formula  from  Section  2.3. 

For  example,  at  the  end  of  a  series  of  iterations  on  the  three-dimensional  flow  field 
through,  say,  blade  row  1,  the  average-passage  source  term  is  calculated  using  a  discrete  version  of 
Equation  2.15  with  i  =  l.  Next,  a  series  of  iterations  is  calculated  for  the  three-dimensional  flow 
through  blade  row  2  utilizing  the  average-passage  source  term  just  calculated  from  the  result  for 
blade  row  1.  Upon  completion  of  these  iterations,  the  average-passage  source  term  for  blade  row  2 
is  computed  using  the  discretized  Equation  2.15  with  i=2.  Blade  row  1  is  then  restarted  with  this 
updated  source  term  and  the  cycle  repeats.  This  periodic  "flipping"  between  blade  rows  using 
updated  average-passage  source  terms  is  what  is  referred  to  as  the  outer  iterations.  In  effect,  at  the 
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4.  Comparison  of  Numerical  and  Experimental  Results 


The  previous  chapters  have  described  the  average-passage  equations  for  a  multiple-blade- 
row  turbomachine  and  the  numerical  technique  used  to  solve  these  equations.  This  description 
includes  a  method  to  model  the  closure  problem  associated  with  a  two-blade-row  machine.  In  this 
chapter,  we  use  the  resulting  computer  code  (ISTAGE)  to  predict  the  flow  through  a  two-blade-row 
pump.  As  a  first  step  towards  code  validation,  this  chapter  presents  comparisons  between  these 
numerical  predictions  and  experimentally  measured  data. 


4.1  High  Reynolds  Number  Pump  Experiment 

Zierke,  Straka,  and  Taylor  [1993]  acquired  a  large  quantity  of  data  within  their  high 
Reynolds  number  pump  (HIREP)  facility.  Figure  4.1  shows  a  computer-generated  graphical  image 
of  the  HIREP  blades  and  the  cylindrical  coordinate  system  used  for  analyzing  the  experimental 
data.  This  isometric  view  shows  that  HIREP  contains  13  inlet  guide  vanes  and  7  rotor  blades  with 
a  significant  circumferential  blade  lean  or  skew.  The  hub  has  a  constant  diameter  of  21  inches, 
while  a  tunnel  liner  creates  a  casing  endwall  with  a  constant  diameter  of  42  inches.  Using  this 
nominal  blade  diameter,  Zierke,  Straka,  and  Taylor  [1993]  found  a  rotor  blade  tip  speed  of 
Utip  =  47.6  ft/sec  while  operating  at  the  design  point  with  260  rpm.  Most  of  the  data  were 
acquired  at  this  design  point.  The  experimental  measurement  techniques  included  rotor  shaft 
torque  and  thrust  cells,  blade  static-pressure  taps,  oil-paint  surface  flow  visualization,  cavitation 
visualization,  laser  light  sheet  visualization,  slow-  and  fast-response  pressure  probes,  and  a  two- 
component,  laser  Doppler  velocimeter  (LDV).  Five-hole  pressure  probe  measurements  of  the 
inflow,  37%  chord  axially  upstream  of  the  inlet  guide  vanes,  showed  a  nearly  uniform,  axial  flow 
with  a  nominal  reference  velocity  of  35  ft/sec.  However,  the  measurements  indicated  that  the 
casing  endwall  boundary  layer  was  fully  turbulent  with  a  boundary  layer  thickness  of  0.71  inches. 
The  turbulent  boundary  layer  on  the  hub  was  slightly  less  than  0.3  inches. 

Since  this  experiment  involved  the  incompressible  flow  of  water  through  a  two-blade-row 
machine,  the  resulting  database  yields  an  excellent  test  case  for  this  particular  code.  In  addition, 
with  chord  Reynolds  numbers  as  high  as  5,500,000,  the  database  also  gives  a  very  challenging  test 
case. 


4.2  Grid  Generation 

Using  the  algebraic  grid  generation  code  of  Beach  [1990]  and  Beach  and  Hoffman  [1992], 
we  generated  a  computational  grid  for  the  single-stage  HIREP  geometry.  This  grid  generator  was 
written  specifically  for  the  compressible  average-passage  formulation.  Basically,  this  interactive 
grid  generator  creates  stacked,  body-conforming,  H-type  grids  for  multiple-blade-row 
turbomachines. 

A  two-blade-row  average-passage  calculation  (like  HIREP)  requires  a  total  of  three  grids. 
There  are  two  three-dimensional  grids:  one  for  an  inlet  guide  vane  (IGV)  passage  and  one  for  a 
rotor  blade  passage.  Each  of  these  grids  extends  from  an  inlet  plane  upstream  of  the  inlet  guide 
vanes  to  an  exit  plane  downstream  of  the  rotor  blades.  Figures  4.2  and  4.3  shows  slices  of  these 
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grids.  The  three-dimensional  grids  used  in  this  work  contain  141  axial  grid  points  and  45  points  in 
both  the  radial  and  tangential  directions.  Additionally,  a  common  two-dimensional  grid  exists  that 
defines  the  meridional  domain  (that  is,  the  extent  in  the  r,z-plane)  of  the  turbomachine.  It  also 
contains  141  axial  and  45  radial  grid  points.  Figure  4.4  shows  this  two-dimensional  grid  for  the 
HIREP  geometry.  It  is  through  this  plane  that  the  two  three-dimensional  solutions  "communicate." 
Body  forces  and  correlations  computed  from  the  three-dimensional  solution  for  one  blade  row  are 
computed  on  this  plane  and  are  subsequently  introduced  into  the  three-dimensional  domain  of  the 
other  blade  row.  Kirtley,  Beach,  and  Adamczyk  [1990]  describe  the  interpolation  procedure 
necessary  to  accomplish  this  interaction. 

In  this  application,  we  generated  the  three-dimensional  grid  for  the  rotor  as  if  the  blade 
physically  extended  from  the  hub  all  the  way  to  the  casing  with  no  tip  clearance.  In  reality,  the 
HIREP  rotor  included  a  nominal  tip  clearance  of  approximately  0.13  inches.  To  model  the 
clearance  flow,  therefore,  we  chose  a  grid  index  that  most  closely  approximated  the  location  of  the 
physical  blade  tip.  A  periodic  boundary  condition  was  then  used  for  all  eight  cells  between  the  tip 
index  and  the  casing.  In  other  words,  above  the  rotor  blade  tip,  a  periodic  condition  was  imposed 
across  the  thickness  of  the  blade  tip  section.  Note  that  we  chose  the  tip  grid  index  to  correspond  to 
the  physical  blade  tip  and  did  not  include  a  discharge  coefficient  to  model  the  vena  contracta  that 
occurs  when  the  tip  clearance  flow  separates.  The  HIREP  rotor  blades  had  a  rounded  pressure  side 
comer  at  the  blade  tip  to  prevent  this  separation  and,  thus,  improve  gap  cavitation  performance. 

In  reality,  both  the  inlet  guide  vanes  and  the  rotor  blades  have  fillets  at  their  junctures  with 
the  hub  and  casing.  Each  inlet  guide  vane,  therefore,  includes  a  fillet  at  its  root  and  tip,  while 
each  rotor  blade  includes  a  fillet  at  the  root  only  (while  the  tip  includes  the  rounded  pressure  side 
comer).  These  fillets,  however,  are  absent  from  the  computational  grid.  Neither  gridding  the  true 
blade  tip  nor  gridding  fillets  is  technically  beyond  the  capabilities  of  the  grid  generation  code; 
however,  making  the  simplifications  just  described  does  greatly  simplify  the  grid  generation 
process.  Moreover,  neither  geometric  simplification  should  introduce  undue  error  into  the 
computational  results. 


4.3  Solution  Procedure 

The  major  assumption  in  the  simulation  of  HIREP  was  the  state  of  the  incoming  flow. 
Without  prior  knowledge  of  this  inflow,  we  used  the  design  inflow  (from  a  streamline  curvature 
simulation).  The  design  inflow  was  simply  uniform  axial  flow.  Zierke,  Straka,  and  Taylor  [1993] 
subsequently  found  from  measurements  of  the  actual  flow  field  that  a  relatively  thick  boundary 
layer  existed  on  the  casing  endwall  upstream  of  the  inlet  guide  vanes.  It  is  believed  that  the 
absence  of  this  boundary  layer  in  the  simulation  may  be  quite  important  to  some  of  the  results-an 
effect  that  we  will  later  address. 

Actually  obtaining  the  average-passage  solution  for  the  HIREP  geometry  involved  running 
the  computer  code  until  a  steady-state  solution  was  reasonably  approximated.  Running  the  code 
involved  a  series  of  "flips"  between  the  two  blade  rows.  Each  "flip"  for  a  given  blade  row 
consisted  of  the  calculation  of  a  number  of  time  steps  (perhaps  several  hundred)  on  the  three- 
dimensional  grid  for  that  blade  row  using  the  (axisymmetric)  average-passage  body  forces  and 
correlations  from  the  other  blade  row  from  the  previous  "flip."  During  the  "flipping"  of  the 
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solutions,  we  monitored  the  mass  flow.  If  the  mass  flow  through  HIREP  was  too  high  or  too  low, 
we  adjusted  the  exit  hub  static  pressure  (back  pressure)  accordingly.  This  "throttling"  of  the 
simulation  was  a  direct  result  of  the  total  pressure/absolute  flow  angle  inflow  boundary  condition 
discussed  in  the  previous  chapter.  The  simulation  was  considered  converged  when  the  calculated 
mass  flow  equaled  the  design  mass  flow  and  both  blade  rows  had  converged  to  the  same 
circumferentially-averaged  flow  field  (recall  that  this  was  a  necessary  convergence  criterion  for  the 
average-passage  equations).  For  the  HIREP  simulation,  the  calculated  mass  flow  was  about  1% 
lower  than  the  design  mass  flow. 


4.4  Comparisons 

Using  the  computational  grid  generated  for  the  HIREP  geometry  and  the  solution  procedure 
just  described,  we  obtained  a  converged  solution  for  the  HIREP  flow  field.  First,  we  shall  discuss 
the  simulation  of  the  IGV  flow  field-including  comparisons  with  the  experimental  database.  Then, 
we  shall  examine  the  flow  field  in  the  vicinity  of  the  downstream  rotor  blades. 


4.4.1  Inlet  Guide  Vanes 

The  first  step  in  understanding  the  basic  characteristics  of  a  turbomachine  flow  field 
normally  involves  observation  of  the  blade  static-pressure  distribution.  For  the  IGV  static 
pressures,  Figure  4.5  shows  a  comparison  between  numerical  and  experimental  results  at  five 
spanwise  locations.  Note  that,  in  this  figure,  the  pressure  coefficient,  Kp,  is  normalized  by  a 
dynamic  pressure  based  on  the  blade  tip  speed— as  defined  by  Zierke,  Straka,  and  Taylor  [1993]. 
Overall,  the  integrated  blade  loading  agrees  reasonably  well  at  all  five  locations.  As  each  section 
is  examined  individually,  we  observe  some  differences  in  the  predicted  and  measured  distributions 
of  pressure. 

At  10%  span,  a  discrepancy  exists  between  the  prediction  and  experiment  in  the 
distribution  of  pressure  near  the  leading  edge.  At  this  location,  lower  pressures  were  measured 
than  were  predicted  on  both  the  suction  and  pressure  surfaces.  Note  that  a  similar  discrepancy  at 
this  spanwise  location  was  observed  when  Zierke,  Straka,  and  Taylor  [1993]  compared  the 
measurements  to  a  lifting  surface  theory  calculation  of  the  flow  field.  In  the  midspan  regions 
(from  30%  to  70%  span),  the  suction  surface  pressures  compare  quite  well  with  the  measurements. 
However,  at  these  same  locations,  on  the  pressure  surface,  the  predictions  show  a  distinct  favorable 
pressure  gradient  near  the  trailing  edge  that  is  not  present  in  the  measurements.  This  phenomenon 
reverses  itself  at  the  90%  span  location,  where  a  favorable  pressure  gradient  was  measured  and  an 
adverse  pressure  gradient  was  predicted. 

The  contour  plots  of  Figures  4.6  and  4.7  reinforce  these  previous  observations.  Notice  that 
the  contours  of  measured  static  pressures  are  not  extrapolated  beyond  the  locations  of  the  pressure 
taps.  Figure  4.6  shows  good  qualitative  agreement  over  the  entire  suction  surface.  For  the 
pressure  surface,  on  the  other  hand,  Figure  4.7  shows  good  agreement  to  about  75%  or  80% 
chord;  beyond  that  location,  the  discrepancies  mentioned  previously  are  clearly  evident.  The 
disagreement  near  the  trailing  edge  on  the  pressure  surface  is  most  likely  due  to  some  unsteadiness 
in  the  flow  solution  there.  The  trailing  edge  geometry  of  the  inlet  guide  vane  is  essentially  an 
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asymmetric  chisel,  with  the  bevel  of  the  chisel  lying  on  the  suction  side  of  the  blade  and  the  point 
of  the  chisel  being  somewhat  rounded  off.  As  such,  a  small  separation  region  exists  in  the  bevel 
for  the  entire  span  of  the  blade.  Depending  on  the  behavior  of  the  wall  function  and  turbulence 
model  in  this  region,  the  predictions  may  quite  possibly  give  an  improperly  located  rear  stagnation 
point.  The  lack  of  a  sharp  trailing  edge  exacerbates  this  potential  problem. 

In  the  IGV  wake  region,  Figure  4.8  shows  the  velocity  distributions  at  twelve  spanwise 
locations  located  49.7%  chord  axially  downstream  of  the  IGV  trailing  edge.  The  figure  compares 
the  axial,  radial,  and  tangential  velocity  components  from  the  prediction  with  those  measured  using 
a  rake  of  five-hole  probes.  In  order  to  make  this  comparison,  we  needed  to  perform  a  three- 
dimensional  interpolation  of  the  numerical  results,  since  the  grid  lines  did  not  lie  in  axial  planes. 

At  4.8%  span  and,  to  a  lesser  extent,  at  9.5%  span,  the  axial  velocity  comparison  shows  that  the 
predicted  wake  depth  and  width  are  smaller  than  the  measurements.  Examination  of  the  predicted 
flow  field  showed  the  comer  separation  on  the  suction  surface  to  be  smaller  than  the  separation  that 
was  measured.  In  the  predicted  flow  field,  this  difference  manifested  itself  downstream  as  a 
shallower,  narrower  wake.  Also,  to  conserve  mass,  the  peak  in  the  predicted  axial  velocity  on  the 
pressure  side  of  the  wake  is  greater  than  the  measured  peak. 

At  larger  spanwise  locations.  Figure  4.8  shows  that  the  numerical  predictions  of  axial 
velocity  give  deeper  wakes  than  the  five-hole  probe  measurements.  This  result  is  surprising  since 
the  grid  points  at  this  axial  location  were  more  sparse  than  near  the  blade  itself-which  should  yield 
larger  numerical  dissipation  and  a  faster  decay  of  the  viscous  wakes.  From  19.0%  span  to  76.2% 
span,  the  axial  velocity  comparisons  show  another  surprise:  A  peak  in  the  predicted  axial  velocity 
occurs  on  the  suction  side  of  the  wake--a  peak  that  does  not  occur  in  the  measurements.  This  peak 
appears  outside  of  the  viscous  wake;  and  yet,  this  flow  field  location  is  too  far  downstream  of  these 
lightly-loaded  inlet  guide  vanes  to  feel  any  potential  flow  effects.  The  deeper,  predicted  wakes 
continue  until  the  95.2%  span  location,  where  the  contamination  of  the  measured  wakes  from 
vortices  emanating  from  axially-running  seams  between  tunnel  liner  sections  made  comparisons 
difficult. 

The  circumferential  variation  of  radial  velocity  in  Figure  4.8  occurs  because  the  spanwise 
variation  in  circulation  on  the  inlet  guide  vanes  gives  rise  to  a  trailing  vortex  sheet.  Avoiding  the 
comparisons  at  4.8%,  90.5%,  and  95.2%  span  where  calibration  problems  produced  offsets  in  the 
radial  velocity  measurements,  we  see  that  the  numerical  predictions  compare  very  well  with  the 
five-hole  probe  measurements  of  the  radial  velocity  variations.  Since  the  trailing  vortex  sheet  is 
radially  skewed,  the  sheet  also  leads  to  a  circumferential  variation  in  tangential  velocity.  The  swirl 
produced  by  the  inlet  guide  vanes  also  means  that  a  small  portion  of  the  viscous  wake  profile  will 
have  a  component  of  tangential  velocity.  Figure  4.8  shows  the  resulting  variation  in  tangential 
velocity.  Despite  the  measured  blade-to-blade  differences,  Figure  4.8  shows  good  comparisons 
between  predicted  and  measured  tangential  velocities.  Near  the  hub,  the  level  of  tangential  velocity 
shows  that  the  predictions  give  a  little  more  swirl  than  what  was  measured. 

The  circumferential  variations  of  radial  and  tangential  velocity  can  also  be  shown  as 
secondary  velocity  vectors.  Figure  4.9  presents  the  flow  within  a  plane  normal  to  the  axial 
direction  with  the  local  circumferential-mean  velocity  (found  from  area  averaging)  subtracted  from 
both  the  predicted  and  measured  velocities.  The  resulting  secondary  velocity  vectors  follow  closely 
to  those  defined  by  Smith  [1955],  Both  the  predictions  and  measurements  show  that  the  trailing 
vortex  sheet  has  deformed  and  rolled-up  into  two  concentrated  vortex  structures  rotating  in  opposite 
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directions.  The  predicted  locations  of  these  two  structures  lies  very  close  to  the  measured 
locations,  although  the  measured  vectors  appear  to  be  somewhat  larger.  These  two  vortex 
structures  induce  a  swirling  motion  in  the  potential  portion  of  the  passage  flow.  Near  the  hub,  the 
predicted  location  of  this  clockwise  (looking  upstream)  motion  matches  the  measured  location;  but 
near  the  casing,  the  predicted  location  of  the  counter-clockwise  motion  appears  closer  to  the 
pressure  side  of  the  wake  than  the  measured  location. 

At  this  same  axial  location  within  HIREP,  Figure  4.10  presents  predicted  and  measured 
contours  of  total-pressure  coefficient.  As  with  the  secondary  velocity  vector  plots,  the  predicted 
total-pressure  contours  are  smoother  than  the  measured  plots  since  the  numerical  simulation 
contained  more  grid  points  than  the  experiment  contained  five-hole  probe  locations.  Both  the 
predictions  and  measurements  showed  a  slightly  radially-skewed  wake  with  similar  levels  of  total- 
pressure  coefficient.  Without  the  casing  boundary  layer  used  in  the  simulation  inflow,  the 
predictions  do  not  show  the  correct  total-pressure  loss  near  the  casing  endwall.  This  effect  is  even 
more  evident  after  performing  a  circumferential  average,  as  shown  in  Figure  4.11.  This  figure 
shows  higher  levels  of  predicted  total-pressure  coefficient  near  the  casing  endwall,  but  lower 
predicted  levels  elsewhere.  The  loss  in  total-pressure  coefficient  through  most  of  the  span  (from  an 
inlet  value  of  unity)  should  be  minimal  through  these  lightly-loaded  vanes.  Note  that  the 
measurements  using  the  rake  of  five-hole  probes  shows  that  the  values  of  total-pressure  coefficient 
remain  close  to  unity,  while  the  measurements  using  the  five-hole  probe  radial  surveys  actually 
show  values  greater  than  unity.  Zierke,  Straka,  and  Taylor  [1993]  explain  that  the  experimental 
uncertainty  is  much  greater  for  the  radial  surveys  than  for  the  rake  surveys.  Figure  4.11  also 
shows  a  larger  predicted  gradient  in  the  static-pressure  profile  than  shown  in  the  measurements  of 
static  pressure— a  quantity  difficult  to  measure  away  from  solid  surfaces. 

Finally,  Figure  4.12  shows  comparisons  of  predicted  and  measured  values  of  the 
circumferentially-averaged  velocity  components.  Overall,  these  comparisons  are  quite  good. 
However,  Figure  4.12  clearly  shows  the  absence  of  the  correct  inflow  boundary  layer  on  the  casing 
endwall  during  the  numerical  simulation.  Consequently,  in  order  to  preserve  the  same  mass  flow, 
the  numerical  simulation  results  in  smaller  values  of  circumferentially-averaged  axial  velocities 
away  from  the  endwalls. 


4.4.2  Rotor  Blades 

Recall  that  the  solution  procedure  used  in  the  numerical  simulation  of  HIREP  consisted  of  a 
series  of  "flips"  through  which  the  three-dimensional  solution  of  the  IGV  passage  "communicated" 
with  the  three-dimensional  solution  of  the  rotor  blade  passage.  This  "communication"  occurred  on 
a  two-dimensional  grid  that  contained  updated  axisymmetric  body  forces  and  correlations  from  the 
latest  "flip"  from  one  blade  row  solution  to  the  other.  Through  this  procedure,  the  momentum  and 
vorticity  field  from  the  IGV  flow  field  solution  described  in  the  previous  section  is  correctly 
transferred  to  the  rotor  blade  flow  field  solution.  And  even  though  the  momentum  and  vorticity 
field  from  the  IGV  flow  field  appears  "smeared"  to  the  rotor  blades,  the  transfer  of  momentum  and 
vorticity  that  does  take  place  should  adequately  model  the  effect  of  the  inlet  guide  vanes  on  the 
time-average  flow  through  the  rotor  blades.  The  procedure  also  adequately  models  the  potential 
flow  effect  of  the  rotor  blades  on  the  time-average  flow  through  the  inlet  guide  vanes. 
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Similar  to  our  comparisons  for  the  IGV  flow  field,  we  will  initially  compare  the  predicted 
and  measured  rotor  blade  static-pressure  distributions.  Figure  4.13  shows  that  these  comparisons 
for  the  rotor  blades  are  even  better  than  the  comparisons  for  the  inlet  guide  vanes  in  Figure  4.5— 
the  rotor  blade  flow  field  showed  less  unsteadiness  in  the  predictions  near  the  trailing  edge  than  did 
the  IGV  flow  field.  Again  using  the  pressure  coefficient,  Kp,  Figure  4.13  shows  that  both  the 
blade  loading  and  the  local  streamwise  pressure  gradients  are  well  predicted  at  all  spanwise 
locations.  The  contour  plots  of  Figures  4.14  and  4.15  also  show  very  good  agreement,  especially 
when  only  comparing  the  predicted  static-pressure  contours  within  the  regions  where  the  measured 
contours  are  bounded  by  the  locations  of  the  pressure  taps.  The  suction  surface  static  pressures  of 
Figure  4.14  show  an  adverse  pressure  gradient  over  the  last  half  of  the  blade,  with  the  predictions 
showing  some  three-dimensional  effects  in  the  comer  region  were  the  trailing  edge  meets  the  hub 
endwall-a  region  where  no  pressure  taps  existed.  Figure  4.15  shows  that  the  predictions  also  give 
some  three-dimensional  effects  where  the  pressure  surface  trailing  edge  meets  the  hub  endwall. 

Numerically  integrating  the  predicted  three-dimensional  static-pressure  distribution  allows 
us  to  calculate  net  thrust  and  torque  coefficients.  Figure  4.16  shows  these  calculations  compared  to 
the  measured  and  design  values.  At  the  design  volumetric  flow  coefficient,  the  predicted  torque 
coefficient  matches  both  the  measured  and  design  values  very  well,  with  the  predicted  value  being 
about  0.7%  higher  than  the  design  value.  (The  measured  torque  coefficient  was  2%  lower  than  the 
design  value.)  Comparisons  of  the  thrust  coefficient  show  that  the  predicted  value  is  8.5%  higher 
than  the  design  value.  While  the  actual  thrust  coefficient  is  probably  higher  than  the  design  value, 
the  actual  coefficient  is  most  likely  less  than  the  measured  value  (which  was  16%  higher  than  the 
design  value).  Zierke,  Straka,  and  Taylor  [1993]  describe  the  difficulties  that  they  had  in 
performing  this  measurement-difficulties  that  led  to  greater  measurement  uncertainties.  Therefore, 
the  simulation  appears  to  give  very  reasonable  predictions  of  rotor  thrust  since  the  predicted  value 
in  Figure  4.16  lies  between  the  measured  and  design  values. 

The  significance  of  these  types  of  numerical  simulations  lies  in  the  ability  to  give  three- 
dimensional,  viscous  flow  information  to  the  designer.  One  very  prominent  characteristic  of  this 
flow  field  prediction  is  the  existence  and  location  of  regions  of  boundary  layer  separation.  The 
next  series  of  figures  shows  how  well  this  numerical  simulation  predicts  two  types  of  three- 
dimensional  separation. 

First,  Figure  4.17  presents  the  predicted  particle  paths  near  the  rotor  hub  endwall. 
Restricted  to  one  grid  point  away  from  the  hub,  these  particle  paths  follow  the  relative  flow 
through  the  rotor  blade  passage  and  can  be  used  to  simulate  limiting  streamlines.  These  limiting 
streamlines  can  be  compared  to  the  skin-friction  lines  found  experimentally  using  an  oil-paint 
technique  and  shown  schematically  in  Figure  4.18.  Recall  that  limiting  streamlines  and  skin- 
friction  lines  are  indistinguishable,  except  for  some  differences  that  occur  in  the  vicinity  of  lines  of 
separation.  Both  Figure  4.17  and  Figure  4.18  show  one  type  of  three-dimensional  separation.  Just 
upstream  of  the  rotor  blade  leading  edge,  the  hub  endwall  boundary  layer  separates  at  a  saddle 
point,  with  a  separation  line  passing  through  the  saddle  point  and  forming  a  C-pattem  around  the 
blade  Although  it  is  difficult  to  identify  the  exact  location  of  the  predicted  saddle  point  from 
Figure  4.17,  it  appears  that  the  predicted  saddle  point  is  a  little  closer  to  the  leading  edge  than  the 
measured  saddle  point  in  Figure  4.18.  The  presence  of  the  fillet  (and  possibly  a  thicker  boundary 
layer)  in  the  actual  flow  field  will  cause  the  saddle  point  to  lie  further  upstream  of  the  predicted 
location.  Nevertheless,  the  numerical  simulation  does  show  the  existence  of  the  saddle  point,  as 
well  as  the  measured  movement  of  the  pressure  side  leg  of  the  separation  line  into  the  midpassage 
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region  and  the  measured  movement  of  the  suction  side  leg  of  the  separation  line  up  onto  the  suction 
surface  as  a  result  of  some  local  flow  overturning.  The  prediction  of  this  type  of  three-dimensional 
separation  is  important  since  it  corresponds  to  the  development  of  horseshoe  vortices. 

The  measurements  of  Figure  4.18  also  shows  the  existence  of  a  saddle  point  near  the 
suction  surface  trailing  edge,  as  well  as  a  corresponding  line  of  secondary  separation.  Clearly,  the 
predicted  particle  paths  in  Figure  4.17  also  show  this  secondary  separation  line  and  the  probable 
existence  of  the  saddle  point.  As  described  by  Zierke,  Straka,  and  Taylor  [1993],  this  secondary 
separation  is  related  to  the  three-dimensional  comer  separation  that  occurs  where  the  suction 
surface  trailing  edge  meets  the  hub  endwall.  Figure  4.19  shows  the  predicted  particle  paths  of  this 
comer  separation  on  the  suction  surface,  while  Figure  4.20  shows  a  schematic  of  the  skin-friction 
lines  found  from  the  oil-paint  technique.  The  predicted  and  measured  patterns  associated  with  this 
corner  separation  look  remarkably  similar,  although  the  predictions  do  not  show  the  existence  of 
the  spiral  node  as  clearly  as  the  measurements.  Both  the  predictions  and  the  measurements  also 
show  a  large  amount  of  outward  radial  migration  of  the  flow.  However,  the  measurements  do 
show  that  the  separation  line  that  extends  out  to  the  blade  tip  lies  further  upstream  of  the  trailing 
edge  than  predicted  by  the  numerical  simulation. 

Next,  we  shall  compare  the  predictions  with  the  measurements  acquired  on  an  axial  plane 
32.2%  chord  axially  downstream  of  the  rotor  tip  trailing  edge.  As  with  the  flow  downstream  of 
the  inlet  guide  vanes,  we  needed  to  perform  a  three-dimensional  interpolation  of  the  numerical 
results,  since  the  grid  lines  did  not  lie  in  axial  planes.  The  predicted  and  measured  axial  velocity 
contours  of  Figure  4.21  both  show  the  existence  of  the  skewed  rotor  blade  wakes.  The  contour 
levels  within  these  wakes  are  quite  similar,  with  the  nonuniform  nature  of  the  measured  wakes 
arising  from  computing  contours  using  LDV  data  from  coarsely-spaced,  discrete  measurement 
locations.  As  a  matter  of  fact,  all  of  the  predicted  and  measured  contours  compare  well  in  this 
plane,  except  for  one  very  important  region:  the  region  near  the  rotor  tip  leakage  vortex.  The 
measured  tip  leakage  vortex  appears  further  from  the  casing  endwall  than  the  predicted  vortex, 
with  the  potential  flow  effects  from  the  measured  vortex  extending  much  further  into  the  flow  field. 
Comparisons  of  the  predicted  and  measured  tangential  velocity  contours  in  Figure  4.22 
quantitatively  give  the  same  conclusions  as  the  comparisons  of  the  axial  velocities 

Similar  to  the  numerical  results  downstream  of  the  inlet  guide  vanes,  we  can  plot  the 
velocities  downstream  of  the  rotor  blades  using  secondary  velocity  vectors,  with  the  local 
circumferential-mean  velocity  (found  from  area  averaging)  subtracted  from  the  predicted  velocities. 
Figure  4.23  presents  this  vector  plot.  No  comparisons  can  be  made  of  these  predicted  vectors  with 
measured  vectors  since  Zierke,  Straka,  and  Taylor  [1973]  could  not  measure  radial  velocities  with 
their  two-component  LDV  system.  Similar  to  the  IGV  trailing  vortex  sheet,  the  rotor  blade  trailing 
vortex  sheet  has  deformed  and  rolled-up  into  two  concentrated  vortex  structures  rotating  in  opposite 
directions.  While  much  of  the  swirling  motion  between  vortex  sheets  has  been  induced  by  the 
sheets  themselves,  one  can  clearly  observe  the  existence  of  the  rotor  tip  leakage  vortex,  which 
again  lies  very  close  to  the  casing  endwall. 

To  more  closely  examine  the  position  of  the  tip  leakage  vortex,  we  interpolated  the 
numerical  results  at  several  axial  planes  downstream  of  the  rotor  tip  trailing  edge.  In  order  to 
identify  the  location  of  the  vortex  core,  we  plotted  contours  of  static  pressure  on  these  interpolated 
planes,  with  the  clearly  marked  position  of  the  minimum  static  pressure  giving  the  location  of  the 
core.  To  summarize,  Figure  4.24  shows  the  spanwise  location  of  the  vortex  core  compared  to 
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measured  values  found  from  both  LDV  data  and  laser  light  sheet  visualization.  This  figure  shows 
that  the  predicted  core  lies  2-3%  span  closer  to  the  casing  endwall,  with  the  radial  locations  for 
both  the  predicted  and  measured  cores  changing  very  little  downstream  of  the  trailing  edge. 

Since  the  predictions  did  not  include  an  inlet  boundary  layer,  the  numerical  simulation 
computed  a  smaller  casing  endwall  boundary  layer  upstream  of  the  rotor  blades  than  shown 
experimentally  (see  Figure  4.12).  If  the  numerical  simulation  had  computed  a  thicker  boundary 
layer,  the  incidence  angle  at  the  rotor  blade  tip  would  have  been  larger  and,  subsequently,  the 
strength  of  the  rotor  tip  leakage  vortex  would  increase.  Using  the  image  vortex  system  that  Chen, 
Greitzer,  Tan,  and  Marble  [1991]  used  within  a  rotor  blade  passage,  a  stronger  vortex  will  change 
the  strength  of  the  image  vortices.  Finally,  the  stronger  image  vortices  should  induce  a  motion  on 
the  tip  leakage  vortex  that  would  cause  the  position  of  the  tip  leakage  vortex  to  move  further  away 
from  the  casing  endwall  (and  slightly  further  away  from  the  rotor  blade  suction  surface).  Later, 
we  recomputed  the  HIREP  flow  field  using  the  measured  inflow.  Unfortunately,  the  revised 
numerical  simulation  did  not  predict  a  significant  movement  of  the  tip  leakage  vortex  away  from 
the  casing  endwall  and,  thus,  the  position  of  the  vortex  still  disagreed  with  the  measured  position. 

Despite  the  discrepancies  between  the  spanwise  positions  of  the  tip  leakage  vortex, 

Figure  4.25  shows  that  the  positions  of  the  predicted  and  measured  vortex  core  agree  very  well  in 
the  blade-to-blade  plane.  To  better  visualize  the  tip  leakage  vortex  in  the  numerical  simulation,  we 
observed  the  paths  of  particles  that  were  released  at  the  edge  of  the  clearance  region  and  rolled-up 
into  a  vortex.  Observation  of  the  paths  of  particles  released  near  the  trailing-edge  tip  reveals 
another  interesting  phenomenon.  These  trailing-edge  particle  paths  seem  to  roll-up  into  an 
independent  vortex,  that  rotates  with  the  same  sense  as  the  tip  leakage  vortex,  but  originates  closer 
to  the  casing  endwall.  Zierke,  Straka,  and  Taylor  [1993]  identified  this  second  vortex  as  a  trailing- 
edge  separation  vortex,  that  originates  when  the  blade  boundary  layer  separates  near  the  trailing 
edge  in  the  presence  of  radial  outward  flow.  In  a  very  similar  experiment,  Farrell  [1989]  was  able 
to  lower  his  tunnel  pressure  to  a  low  enough  level  for  cavitation  to  occur  in  both  vortices. 

Figure  4.27  shows  his  photograph  of  this  phenomenon-an  observation  very  similar  to  the  predicted 
particle  paths  of  Figure  4.26.  Further,  cavitation  visualization  performed  by  Farrell  [1989]  showed 
that  these  two  vortices  will  eventually  roll-up  into  a  single  vortex  as  they  propagate  downstream  (as 
also  indicated  by  the  predictions  in  Figure  4.2b). 

For  a  more  detailed  comparison  between  the  numerical  simulation  and  the  experiment, 
Figure  4.28  presents  15  spanwise  plots  of  the  circumferential  variations  of  axial  and  tangential 
velocity  at  the  axial  plane  32.2%  chord  axially  downstream  of  the  rotor  tip  trailing  edge.  At  2.9% 
span,  the  predictions  match  the  LDV  data  fairly  well,  although  the  experiment  indicates  a  flow 
structure  on  either  side  of  the  wake  (possibly  the  horseshoe  vortex)  that  is  not  predicted.  At  larger 
spanwise  locations,  through  57.7%  span,  the  predictions  agree  fairly  well  with  the  LDV  data, 
although  a  few  differences  are  worth  noting.  At  all  these  spanwise  locations,  the  predicted  viscous 
wakes  are  deeper  than  the  measured  wakes.  Also,  by  observing  the  centers  of  the  wakes,  the 
predicted  wakes  appear  to  be  slightly  more  skewed  than  the  measured  wakes.  Finally,  from  29.4% 
span  to  57.7%  span,  the  overall  level  of  tangential  velocity  is  more  negative  for  the  predictions, 
indicating  the  numerical  simulation  allows  for  more  flow  turning  than  the  experiment. 

Closer  to  the  casing  endwall,  Figure  4.28  shows  poorer  agreement  between  the  predictions 
and  the  measurements.  First,  as  clearly  shown  from  81.3%  span  to  90.7%  span,  the  predicted 
wakes  are  shallower  than  the  measured  wakes-the  opposite  trend  from  what  we  observed  at 
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smaller  spanwise  locations.  During  the  HIREP  experiment,  Zierke,  Straka,  and  Taylor  [1993] 
acquired  LDV  data  at  76.2%  span  and  at  several  axial  locations,  in  order  to  investigate  wake 
decay.  Figure  4.29  presents  comparisons  of  the  predicted  wakes  with  these  LDV  data.  All  the 
measured  wakes  are  narrower  than  the  predicted  wakes.  However,  the  measured  wakes  are  deeper 
at  2.5%  chord  downstream  of  the  trailing  edge  and  reach  a  similar  depth  at  16.5%  chord 
downstream— indicating  a  faster  decay  rate  for  the  experimental  wakes.  Recall  that  in  the  numerical 
simulation,  the  eddy  viscosity  downstream  of  the  blades  is  extrapolated  along  axial  grid  lines.  If 
the  axial  grid  lines  do  not  reasonably  follow  the  wake,  then  the  predicted  eddy  viscosity  in  the 
wake  will  be  too  small.  Therefore,  close  to  the  trailing  edge,  one  might  expect  reasonable 
agreement.  However,  as  the  predicted  wake  moves  downstream  and  strays  from  the  region  where 
the  eddy  viscosity  exists,  one  might  expect  the  wake  depth  to  decay  more  like  a  laminar  wake 
rather  than  a  turbulent  wake. 

Outside  of  the  viscous  wakes,  in  an  essentially  potential  flow  core,  Figure  4.29  shows 
significant  differences  between  the  predictions  and  the  measurements.  Figure  4.28  shows  these 
same  differences  from  67.1%  span  to  86.0%  span,  the  measured  axial  velocities  show  a  greater 
hump,  while  the  measured  tangential  velocities  show  a  greater  trough.  These  flow  structures  result 
from  potential  flow  effects  caused  by  the  rotor  tip  leakage  vortex  and,  since  the  predicted  vortex 
lies  closer  to  the  casing  endwall,  the  potential  flow  effect  from  the  predicted  vortex  will  not  extend 
as  far  inboard.  From  90.7%  span  to  99.2%  span,  the  discrepancies  in  the  radial  position  of  the  tip 
leakage  vortex  core  becomes  even  more  apparent. 

The  incorrect  predicted  radial  position  of  the  tip  leakage  vortex  also  impacts  the 
circumferentially-averaged  velocity  profiles  in  Figure  4.30.  Near  the  casing  endwall,  the  deficit  in 
axial  velocity  results  more  from  circumferentially  averaging  through  the  axial  velocity  deficits  in 
the  tip  leakage  vortex  than  from  the  endwall  boundary  layer.  With  the  predicted  vortex  positioned 
nearer  to  the  casing  endwall,  Figure  4.30  shows  that  the  axial  velocity  deficit  lies  closer  to  the 
casing  for  the  predictions  than  for  the  measurements.  This  predicted  deficit  displaces  less  fluid 
away  from  the  casing  endwall  and  results  in  a  smaller  axial  velocity  through  most  of  the  span  in 
order  to  have  the  same  mass  flow.  Zierke,  Straka,  and  Taylor  [1993]  point  out  that  their 
experimental  results  do  not  average  through  any  IGV  wakes,  and  this  also  increases  the 
circumferentially-averaged  measurements  relative  to  the  predictions.  The  tangential  velocity  profile 
compares  fairly  well  near  the  endwalls;  but  again,  the  predictions  exhibit  more  negative  tangential 
velocity  through  most  of  the  span.  The  inlet  guide  vanes  in  HIREP  were  designed  to  place  positive 
tangential  velocity  into  the  flow  and  the  rotor  blades  were  designed  to  take  all  of  the  swirl  out  of 
the  flow  at  the  blade  tip  and  leave  some  negative  tangential  velocity  in  the  flow  at  the  blade  root. 
Finally,  the  predictions  agree  with  the  HIREP  design  in  that  the  predicted  level  of 
circumferentially-averaged  radial  velocity  is  zero  for  the  entire  span. 

Figure  4.31  shows  the  circumferentially-averaged  predictions  of  the  static-  and  total- 
pressure  coefficients.  The  numerical  simulation  predicts  a  constant  static  pressure  over  the  entire 
span,  with  the  static  pressures  near  the  casing  endwall  predicted  to  be  somewhat  smaller  than  the 
pressure  tap  measurements.  Figure  4.31  also  compares  the  predicted  total  pressures  with  those 
measured  with  radial  Kiel  probe  surveys  at  two  circumferential  positions  of  the  inlet  guide  vanes. 

If  these  experimental  surveys  would  have  been  acquired  for  many  IGV  positions  and  then 
circumferentially  averaged,  these  averaged  values  of  total-pressure  coefficient  would  probably  agree 
quite  well  with  the  predicted  values.  Before  we  performed  the  circumferential  averaging,  the  total- 
pressure  coefficients  varied  as  they  appear  in  the  contour  plot  of  Figure  4.32.  The  higher  regions 
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of  loss  appear  as  low  regions  of  total  pressure  in  the  rotor  blade  wakes  and  near  the  endwalls. 
Figure  4.32  agrees  qualitatively  with  that  data  the  Zierke,  Straka,  and  Taylor  [1993]  measured  with 
a  fast-response  total-pressure  probe.  However,  these  measurements  could  only  measure  total- 
pressure  variations  and  also  included  an  IGV  wake  and,  therefore,  we  could  not  make  a  direct 
comparison  with  the  predictions. 
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5.  Summary  and  Conclusions 


A  system  of  equations  is  presented  governing  the  average-passage  flow  through  a  multiple- 
blade-row  turbomachine  operating  in  the  incompressible  flow  regime.  The  governing  equations 
result  from  the  sequential  application  of  three  averaging  operators  to  the  incompressible,  Navier- 
Stokes  equations.  A  procedure  for  closing  this  averaged  system  of  equations  is  also  described.  To 
facilitate  numerical  treatment,  the  governing  equations  are  cast  in  a  preconditioned,  pseudo¬ 
unsteady  form.  It  is  important  that  the  averaging  operators  be  carried  out  before  casting  the 
equations  in  the  pseudo-unsteady  form.  The  inclusion  of  the  preconditioning  matrix  allows  for 
better  control  of  the  solution  transient. 

An  explicit  numerical  treatment  is  outlined  to  approximate  solutions  to  the  governing 
equations.  Basically,  this  procedure  follows  the  fairly  standard  four-stage  Runge-Kutta  time¬ 
stepping  scheme  applied  to  a  cell-centered  finite-volume  discretization  for  hexahedral  cells.  Several 
convergence  acceleration  techniques-including  local  time  stepping,  implicit  residual  smoothing, 
and  multigrid-are  applied  to  this  scheme  in  order  to  improve  its  efficiency.  Also,  boundary 
conditions  appropriate  for  an  average-passage  simulation  are  presented  Finally,  the  practical 
aspects  of  numerically  closing  the  discretized  governing  equations  are  discussed.  Specifically,  a 
two-tier  iteration  procedure  is  described,  where  the  inner  iterations  update  the  three-dimensional 
solutions  for  each  blade  row  and  the  outer  iterations  update  the  average-passage  source  for  each 
blade  row. 

An  experimental  investigation  is  briefly  described  involving  detailed  measurements  made  of 
the  incompressible  flow  field  in  a  high  Reynolds  number  pump.  This  42- inch  diameter  pump 
includes  two  blade  rows  that  operate  within  a  48-inch  diameter  water  tunnel.  A  numerical 
simulation  using  the  incompressible,  average-passage  equations  was  performed  in  order  to  compare 
with  these  flow  field  measurements. 

Comparisons  between  the  computed  and  measured  static  pressures  on  the  inlet  guide  vanes 
show  reasonable  predictions  of  blade  loading.  However,  some  discrepancies  occur  in  the  pressure 
distributions  near  the  trailing  edge  at  midspan.  The  IGV  wake  velocity  profiles  also  show 
reasonable  comparisons  near  the  endwalls,  with  the  largest  discrepancies  (especially  in  axial 
velocity)  occurring  near  midspan.  Residual  unsteadiness  in  the  numerical  solutions  probably 
account  for  the  discrepancies  in  the  pressure  and  velocity  flow  fields  near  the  IGV  trailing  edge  at 
midspan.  This  unsteadiness  is  most  likely  due  to  the  combination  of  the  bluntness  of  the  IGV 
trailing  edge  and  the  poor  performance  of  the  Baldwin-Lomax  eddy  viscosity  model  in  the  wake 
regions  (at  least  for  this  grid).  Comparison  of  the  total-pressure  contours  and  the 
circumferentially-averaged  pressures  and  velocities  show  good  agreement  away  from  the  casing 
endwall,  where  the  calculation  lacked  the  proper  inflow  boundary  layer. 

Comparisons  of  the  computed  and  measured  rotor  blade  static-pressure  distribution  show 
excellent  agreement  at  all  radial  locations.  As  such,  the  integrated  effect  of  the  computed  rotor 
blade  pressures  gave  good  agreement  with  the  measured  torque  and  thrust.  Computed  particle 
pathlines  on  the  rotor  hub  and  suction  surfaces  qualitatively  captured  the  essential  separation 
features  showed  by  the  experimental  skin-friction  lines.  Inadequate  grid  resolution  probably  led  to 
any  discrepancies  in  these  comparisons.  Comparisons  of  the  contours  of  the  axial  and  tangential 
velocities  downstream  of  the  rotor  blades  show  qualitative  agreement  away  from  the  casing 
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endwall.  Near  the  casing,  the  lack  of  the  correct  inflow  boundary  layer  had  some  detrimental 
effect.  However,  even  this  improper  inflow  boundary  condition  cannot  completely  explain  why  the 
calculated  tip  leakage  vortex  lies  closer  to  the  casing  than  the  measured  vortex.  Despite  this 
discrepancy  in  the  radial  location  of  the  tip  leakage  vortex,  the  blade-to-blade  trajectory  appears  to 
agree  quite  well  with  the  measurements.  The  computations  also  appear  to  capture  the  essential 
character  of  the  measured  leakage  flow;  that  is,  a  primary  vortex  exists,  as  well  as  a  weaker 
trailing-edge  separation  vortex. 

Downstream  of  the  rotor  blades,  comparisons  of  wake  velocity  profiles  show  qualitative 
agreement  away  from  the  endwalls.  Very  close  to  the  hub,  inadequate  grid  resolution  and  the 
subsequent  lack  of  a  well-defined  horseshoe  vortex  in  the  calculation  are  probably  the  culprits. 
Beyond  about  70%  span,  the  discrepancy  in  the  radial  location  of  the  tip  leakage  vortex  causes 
differences  in  the  computed  and  measured  velocity  profiles.  At  76.2%  span,  the  measured  wakes 
decay  more  rapidly  than  the  predicted  wakes.  This  is  most  certainly  a  manifestation  of  the  current 
implementation  of  the  Baldwin-Lomax  model  in  the  wake  region.  Finally,  as  with  the  IGV  flow 
field,  comparisons  of  the  total-pressure  contours  and  the  circumferentially-averaged  pressures  and 
velocities  downstream  of  the  rotor  blades  show  good  agreement  away  from  the  casing  endwall, 
where  the  calculation  lacked  the  proper  inflow  boundary  layer. 

The  comparisons  in  this  investigation  represent  an  initial  validation  case  of  the  numerical 
algorithm  for  the  simulation  of  incompressible  flow  within  a  multiple-blade-row  turbomachine,  as 
outlined  previously.  For  the  most  part,  the  comparisons  indicate  that  the  model  does  an  adequate 
job  of  predicting  many  of  the  essential  flow  features  present  in  a  multiple-blade-row  environment. 
For  example,  blade  pressure  distributions  are  predicted  quite  well-as  well  as  the  overall 
performance,  as  indicated  by  circumferentially-averaged  quantities  and  rotor  blade  torque  and 
thrust.  Also,  certain  highly  three-dimensional  effects  like  leakage  vortices  and  comer  separations 
can  be  qualitatively  predicted.  On  the  other  hand,  it  appears  as  though  the  behavior  of  the 
turbulence  model  in  wake  regions  consistently  causes  problems  in  predicting  wake  velocity  profiles 
(at  least  for  this  grid).  Therefore,  the  most  pressing  need  for  the  numerical  simulation  of  the 
incompressible,  average-passage  flow  through  a  multiple-blade-row  turbomachine  appears  to  be  a 
turbulence  model  that  better  represents  the  physics  in  wake  regions. 
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Figure  3.1  Generic  Flow  Domain  with  the  Relevant  Coordinate  System 
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Figure  3.2  Generic  Finite-Volume  Cell 
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Figure  3.3  Directed  Areas  from  a  Generic  Finite-Volume  Cell 
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Figure  3.4  Neutral  Stability  Boundaries 
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Figure  3.5  Types  and  Locations  of  Boundary  Conditions  for  a  Typical  Blade  Passage 
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4.1  Isometric  View  of  the  HIREP  Blades  with  the  Cylindrical  Coordinate  System 
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Slices  of  the  IGV  Three-Dimensional,  Computational  Grid 
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Slices  of  the  Rotor  Blade  Three-Dimensional,  Computational  Grid 
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Two-Dimensional  Computational  Grid  Defining  the  Meridional  Domain  of  HIREP 
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Figure  4.5  IGV  Static-Pressure  Distribution 
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Figure  4.6  IGV  Suction  Surface  Static-Pressure  Contours 
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Figure  4.7  IGV  Pressure  Surface  Static-Pressure  Contours 
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Figure  4.8  Circumferential  Variation  of  Velocity  Components  49.7%  Chord  Axially  Downstream  of  the 
IGV  Trailing  Edge--Five-Hole  Probe  Measurements  and  ISTAGE  Predictions:  (a)  4.8%  Span 
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Figure  4.8  Circumferential  Variation  of  Velocity  Components  49.7%  Chord  Axially  Downstream  of  the 
IGV  Trailing  Edge--Five-Hole  Probe  Measurements  and  ISTAGE  Predictions:  (b)  9.5%  Span 
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Figure  4.8  Circumferential  Variation  of  Velocity  Components  49.7%  Chord  Axially  Downstream  of  the 

IGV  Trailing  Edge--Five-Hole  Probe  Measurements  and  ISTAGE  Predictions:  (c)  19.0%  Span 
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Figure  4.8  Circumferential  Variation  of  Velocity  Components  49.7%  Chord  Axially  Downstream  of  the 

IGV  Trailing  Edge— Five-Hole  Probe  Measurements  and  ISTAGE  Predictions:  (d)  38.1%  Span 
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Figure  4.8  Circumferential  Variation  of  Velocity  Components  49.7%  Chord  Axially  Downstream  of  the 

IGV  Trailing  Edge--Five-Hole  Probe  Measurements  and  ISTAGE  Predictions:  (e)  52.4%  Span 
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Figure  4.8  Circumferential  Variation  of  Velocity  Components  49.7%  Chord  Axially  Downstream  of  the 
IGV  Trailing  Edge--Five-Hole  Probe  Measurements  and  ISTAGE  Predictions:  (f)  61.9%  Span 
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Figure  4.8  Circumferential  Variation  of  Velocity  Components  49.7%  Chord  Axially  Downstream  of  the 

IGV  Trailing  Edge--Five-Hole  Probe  Measurements  and  ISTAGE  Predictions:  (g)  71.4%  Span 
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Figure  4.8  Circumferential  Variation  of  Velocity  Components  49.7%  Chord  Axially  Downstream  of  the 

IGV  Trailing  Edge--Five-Hole  Probe  Measurements  and  ISTAGE  Predictions:  (h)  76.2%  Span 
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Figure  4.8  Circumferential  Variation  of  Velocity  Components  49.7%  Chord  Axially  Downstream  of  the 
IGV  Trailing  Edge--Five-Hole  Probe  Measurements  and  ISTAGE  Predictions:  (i)  81.0%  Span 
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Figure  4.8  Circumferential  Variation  of  Velocity  Components  49.7%  Chord  Axially  Downstream  of  the 
IGV  Trailing  Edge-Five-Hole  Probe  Measurements  and  ISTAGE  Predictions:  (j)  85.7%  Span 
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Figure  4.8  Circumferential  Variation  of  Velocity  Components  49.7%  Chord  Axially  Downstream  of  the 

IGV  Trailing  Edge--Five-Hole  Probe  Measurements  and  ISTAGE  Predictions:  (k)  90.5%  Span 
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Figure  4.8  Circumferential  Variation  of  Velocity  Components  49.7%  Chord  Axially  Downstream  of  the 
IGV  Trailing  Edge-Five-Hole  Probe  Measurements  and  ISTAGE  Predictions:  (1)  95.2%  Span 
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Figure  4.9  Secondary  Velocity  Vectors  49.7%  Chord  Axially  Downstream 
of  the  IGV  Trailing  Edge:  (a)  ISTAGE  Predictions 
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Figure  4.9  Secondary  Velocity  Vectors  49.7%  Chord  Axially  Downstream 
of  the  IGV  Trailing  Edge:  (b)  Five-Hole  Probe  Rake  Surveys 
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Figure  4.10  Total-Pressure  Coefficient  Contours  49.7%  Chord  Axially  Downstream 
of  the  IGV  Trailing  Edge:  (a)  ISTAGE  Predictions 
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Figure  4.10  Total-Pressure  Coefficient  Contours  49.7%  Chord  Axially  Downstream 
of  the  IGV  Trailing  Edge:  (b)  Five-Hole  Probe  Rake  Surveys 
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Figure  4.11  Circumferentially-Averaged  Pressure  Coefficients  49.7% 
Chord  Axially  Downstream  of  the  IGV  Trailing  Edge 
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Figure  4.12  Circumferentially-Averaged  Velocities  49.7% 

Chord  Axially  Downstream  of  the  IGV  Trailing  Edge 
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Figure  4.13  Rotor  Blade  Static-Pressure  Distribution 
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Figure  4. 14  Rotor  Blade  Suction  Surface  Static-Pressure  Contours 
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Figure  4.15  Rotor  Blade  Pressure  Surface  Static-Pressure  Contours 
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Volumetric  Flow  Coeffic ient 
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4.17  Predicted  Particle  Paths,  One  Grid  Point  from  the  Rotor  Blade  Hub 


4.18  Schematic  of  Surface  Flow  Visualization  on  the  Rotor  Blade  Hub  Surface 
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Figure  4.19  Predicted  Particle  Paths,  One  Grid  Point  from  the  Rotor  Blade  Suction  Surface 


4.20  Schematic  of  Surface  Flow  Visualization  on  the  Rotor  Blade  Suction  Surface 
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Figure  4.21  Axial  Velocity  Contours  32.2%  Chord  Axially  Downstream 
of  the  Rotor  Tip  Trailing  Edge:  (a)  ISTAGE  Predictions 
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Figure  4.21  Axial  Velocity  Contours  32.3%  Chord  Axially  Downstream 
of  the  Rotor  Tip  Trailing  Edge:  (b)  LDV  Measurements 
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Figurfe  4.22  Tangential  Velocity  Contours  32.2%  Chord  Axially  Downstream 
of  the  Rotor  Tip  Trailing  Edge:  (a)  ISTAGE  Predictions 


o  o  o  o 
m  o  u?  o 
o  o 

dodo 


If} 

o 

in 

o 

LO 

o 

r- 

T— 

CM 

CM 

d 

d 

d 

d 

d 

iti  p  i 


> 


36 


Figure  4.22  Tangential  Velocity  Contours  32.2%  Chord  Axially  Downstream 
of  the  Rotor  Tip  Trailing  Edge:  (b)  LDV  Measurements 
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Figure  4.23  Predicted  Secondary  Velocity  Vecotrs  32.2%  Chord  Axially  Downstream  of  the  Rotor  Tip  Trailing  Edge 
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Figure  4.24  Spanwise  Position  of  the  Core  of  the  Rotor  Tip  Leakage  Vortex 
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Figure  4.25  Blade-to-Blade  Position  of  the  Core  of  the  Rotor  Tip  Leakage  Vortex 


Figure  4.26  Prediction  Particle  Paths  of  the  Rotor  Tip  Leakage  Vortex  and  the  Trailing  Edge  Separation  Vortex 
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Figure  4.27  Photograph  of  a  Cavitating  Rotor  Tip  Leakage  Vortex  and  a  Cavitating  Trailing  Edge  Separation  Vortex 
(Farrell  [1989]) 
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Figure  4.28  Circumferential  Variation  of  Velocity  Components  32.2%  Chord  Axially  Downstream 

of  the  Rotor  Tip  Trailing  Edge-LDV  Measurements  and  ISTAGE  Predictions:  (a)  2.9%  Span 
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Figure  4.28  Circumferential  Variation  of  Velocity  Components  32.2%  Chord  Axially  Downstream 

of  the  Rotor  Tip  Trailing  Edge--LDV  Measurements  and  ISTAGE  Predictions:  (b)  8.6%  Span 
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Figure  4.28  Circumferential  Variation  of  Velocity  Components  32.2%  Chord  Axially  Downstream 

of  the  Rotor  Tip  Trailing  Edge--LDV  Measurements  and  ISTAGE  Predictions:  (c)  10.5%  Span 
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Figure  4.28  Circumferential  Variation  of  Velocity  Components  32.2%  Chord  Axially  Downstream 

of  the  Rotor  Tip  Trailing  Edge~LDV  Measurements  and  ISTAGE  Predictions:  (d)  19.9%  Span 
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Figure  4.28  Circumferential  Variation  of  Velocity  Components  32.2%  Chord  Axially  Downstream 

of  the  Rotor  Tip  Trailing  Edge--LDV  Measurements  and  ISTAGE  Predictions:  (e)  29.4%  Span 
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Figure  4.28  Circumferential  Variation  of  Velocity  Components  32.2%  Chord  Axially  Downstream 

of  the  Rotor  Tip  Trailing  Edge-LDV  Measurements  and  ISTAGE  Predictions:  (f)  38.8%  Span 
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Figure  4.28  Circumferential  Variation  of  Velocity  Components  32.2%  Chord  Axially  Downstream 

of  the  Rotor  Tip  Trailing  Edge--LDV  Measurements  and  ISTAGE  Predictions:  (g)  48.2%  Span 
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Figure  4.28  Circumferential  Variation  of  Velocity  Components  32.2 %  Chord  Axially  Downstream 

of  the  Rotor  Tip  Trailing  Edge--LDV  Measurements  and  ISTAGE  Predictions:  (h)  57.7%  Span 
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Figure  4.28  Circumferential  Variation  of  Velocity  Components  32.2%  Chord  Axially  Downstream 

of  the  Rotor  Tip  Trailing  Edge-LDV  Measurements  and  ISTAGE  Predictions:  (i)  67.1%  Span 
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Figure  4.28  Circumferential  Variation  of  Velocity  Components  32.2%  Chord  Axially  Downstream 

of  the  Rotor  Tip  Trailing  Edge— LDV  Measurements  and  ISTAGE  Predictions:  (j)  81.3%  Span 
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Figure  4.28  Circumferential  Variation  of  Velocity  Components  32.2%  Chord  Axially  Downstream 

of  the  Rotor  Tip  Trailing  Edge~LDV  Measurements  and  ISTAGE  Predictions:  (k)  86.0%  Span 
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Figure  4.28  Circumferential  Variation  of  Velocity  Components  32.2%  Chord  Axially  Downstream 

of  the  Rotor  Tip  Trailing  Edge—LDV  Measurements  and  ISTAGE  Predictions:  (1)  90.7%  Span 
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Figure  4.28  Circumferential  Variation  of  Velocity  Components  32.2%  Chord  Axially  Downstream 

of  the  Rotor  Tip  Trailing  Edge--LDV  Measurements  and  ISTAGE  Predictions:  (m)  94.5%  Span 
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Figure  4.28  Circumferential  Variation  of  Velocity  Components  32.2%  Chord  Axially  Downstream 

of  the  Rotor  Tip  Trailing  Edge-LDV  Measurements  and  1ST AGE  Predictions:  (n)  98.3%  Span 
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Figure  4.28  Circumferential  Variation  of  Velocity  Components  32.2%  Chord  Axially  Downstream 

of  the  Rotor  Tip  Trailing  Edge-LDV  Measurements  and  ISTAGE  Predictions:  (o)  99.2%  Span 
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Figure  4.29  Passage-Averaged  Rotor  Blade  Wakes  at  76.2%  Span:  (a)  2.5%,  4.8%,  7.2%,  and  11.8% 
Chord  Axially  Downstream  of  the  Rotor  Blade  Trailing  Edge 
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Figure  4.29  Passage-Averaged  Rotor  Blade  Wakes  at  76.2%  Span:  (b)  16.5%,  21.1%,  and  24.6% 
Chord  Axially  Downstream  of  the  Rotor  Blade  Trailing  Edge 
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Figure  4.30  Circumferentially-Averaged  Velocities  32.2%  Chord  Axially  Downstream  of  the  Rotor  Tip  Trailing  Edge 
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Figure  4.31  Circumferentially-Averaged  Pressure  Coefficients  32.2%  Chord  Axially  Downstream  of  the  Rotor  Tip  Trailing  Edge 
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Figure  4.32  Predicted  Contours  of  Total-Pressure  Coefficient  32.2%  Chord  Axially  Downstream  of  the  Rotor  Tip  Trailing  Edge 
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